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ABSTRACT 


With the recent advances in high-speed digital circuits, modeling of interconnects 
and associated discontinuities has gained a considerable interest over the last decade 
although the theoretical bases for analyzing these structures were well-established as 
early as the 1960s. Ongoing research at the present time is focused on devising methods 
which can be applied to more general geometries than the ones considered in earlier days 
and, at the same time, improving the computational efficiency and accuracy of these 
methods. 

In this thesis, numerically efficient methods to compute the transmission line 
parameters of a multiconductor system and the equivalent capacitances of various strip 
discontinuities are presented based on the quasi-static approximation. The presented 
techniques are applicable to conductors embedded in an arbitrary number of dielectric 
layers with two possible locations of ground planes at the top and bottom of the dielectric 
layers. The cross-sections of conductors can be arbitrary as long as they can be described 
with polygons. 

An integral equation approach in conjunction with the collocation method is used 
in the presented methods. A closed-form Green’s function is derived based on weighted 
real images thus avoiding nested infinite summations in the exact Green’s function; 
therefore, this closed-form Green’s function is numerically more efficient than the exact 
Green’s function. All elements associated with the moment matrix are computed using 
the closed-form formulas. Various numerical examples are considered to verify the 
presented methods, and a comparison of the computed results with other published results 
showed good agreement. 
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CHAPTER 1 
INTRODUCTION 

1.1 Background 

In recent years, the numerical modeling and simulation of interconnections and their 
discontinuities in digital integrated circuits have gained significant interest due to the 
modem development of VLSI technology. As the complexity, density, and speed of the 
integrated circuits continue to increase, signal delay and rise times arc increasingly limited 
by interconnections rather than device speeds, and accurate estimations of the signal delay 
and distortion due to interconnections become crucial at virtually every level of circuit 
integration. 

To accurately characterize signal delay distortion and crosstalk noise due to 
interconnection lines, interconnects must be modeled as multiconductor transmission lines 
instead of conventional lumped circuit elements, and associated discontinuities, such as 
crossovers, bends, junctions, and vias, must also be accurately modeled. Although a 
substantial amount of work has been performed over the last three decades to characterize 
interconnections and their discontinuities in the electromagnetic community [l]-[3], most of 
these theoretical studies resulted in methods which involve high computational cost and, 
hence, are not suitable for the real-time design of CAD tools. 

To overcome this difficulty associated with the theoretical analysis, a model-based 
interconnect capacitance extraction tool is studied in the circuit community [4]-[6]. In the 
model-based approach, analytical or table-look-up models are fitted to the data generated by 
numerical simulation of EM-based techniques or experimental measurements. Although 
this approach may reduce the time to compute parameters associated with interconnects, it 
requires an impractical number of models, which limits its practical usage. Fortunately, 
even when the layout of a circuit is very complex, the number of distinct interconnections 
and their discontinuities is often very limited; furthermore, an accurate characterization of 
interconnects is required only for the critical components (path) in the circuit. Thus, if a 
method based on the electromagnetic analysis is sufficiently fast, it may be incorporated 
into a layout CAD tool. 

This thesis focuses on the discussion of computationally efficient methods for 
interconnection modeling. In particular, this thesis presents methods based on the quasi- 
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static approximation to compute the transmission parameters of a multiconductor inter- 
connection line and the equivalent capacitances of interconnection discontinuities embedded 
in a multilayered dielectric medium. 

1.2 The Quasi-Static Approximation 

A multiconductor transmission line in a multilayered dielectric medium does not 
support TEM modes due to the inhomogeneity of a dielectric medium [7], and full-wave 
analysis must be considered to accurately characterize hybrid modes in the transmission 
line. However, when the transverse components of the electric and magnetic fields are 
predominant over the longitudinal components, the fundamental hybrid mode becomes a 
quasi-TEM mode, in which TEM properties dominate the hybrid modes, and lines which 
support a quasi-TEM mode are called quasi-TEM lines (similarly, lines supporting a TEM 
mode are called TEM lines.) The valid range of a quasi-TEM mode is often determined 
using the dimensional analysis on the Maxwell equations [8], [9]. For most quasi-TEM 
lines, a quasi-TEM mode is valid up to several gigahertz; particularly, it is valid up to the 
cutoff frequency of the next higher (hybrid) mode. 

Since the electric field distribution of a TEM mode is identical to that of the static 
case [ch. 3, 10] and static analysis is simpler and computationally less intensive than full- 
wave analysis, quasi-TEM lines are often analyzed using a static analysis, which is then 
called the quasi-static approximation. All methods presented in this thesis are based on this 
quasi-static approximation. It should be noted that although not all transmission lines are 
quasi-TEM lines, interconnections encountered in digital integrated circuits belong to quasi- 
TEM lines 

1.3 Electrostatic Solution Techniques 

Under the quasi-static approximation, the analysis of interconnects and 
discontinuities is performed by solving electrostatic and magnetostatic problems. As will 
be discussed in Chapter 3, for two-dimensional problems, for example, solving for the 
transmission parameters of interconnects, there exists an analogy between electrostatic and 
magnetostatic problems; therefore, the solutions of two-dimensional magnetostatic 
problems can be obtained by solving the equivalent electrostatic problems. Furthermore, 
since this thesis focuses on the modeling of only the capacitive nature of discontinuities for 
three-dimensional problems, electrostatic solution techniques are sufficient for analyzing 
the problems considered in this thesis. 
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Electrostatic problems are governed by Laplace’s equation or Poisson’s equation 
with appropriate boundary conditions. Various methods have been employed to obtain the 
solution in two-dimensional space [11 ]-[19]. Two commonly used techniques for both 
2-D and 3-D are integral equation methods [20], [21] and the domain methods, such as the 
finite difference method (FD) [22] and finite element method (FEM) [23], [24]. 

In the domain methods, the unknown potential distribution is solved to compute the 
charge distribution over an entire domain by either directly approximating the differential 
equation with the finite difference equation (FD) or using the equivalent variational 
expression in conjunction with the method of subareas (FEM). The major disadvantage of 
the domain methods is that the unknown potential distribution to be sought is over the 
entire geometry considered, including the dielectric region; hence, it may be 
computationally inefficient for the open geometry case even with the use of absorbing 
boundary conditions to truncate the open geometry. The computational inefficiency of the 
finite difference method is improved by employing the method of line (MoL) [25], [26]. In 
MoL, all but one of the independent variables of Laplace’s equation are discretized to obtain 
a system of ordinary linear differential equations. These ordinary differential equations are 
then decoupled using the orthogonal transformation matrix and solved analytically. 
Although MoL is computationally very efficient for two-dimensional problems, it is still 
burdensome for three-dimensional problems. Moreover, this method is only applicable to 
infinitely thin conductors. 

On the other hand, the conventional integral equation approach first obtains the 
Green’s function for a layered medium using the image theory, which consists of rather 
slowly converging infinite series. Then, an integral equation is formulated using this 
Green’s function as its kernel and is solved by employing the method of moments (MoM) 
to determine the unknown charge density on the conductor surfaces. Since unknowns in 
this approach only lie on the surface of conductors, methods based on an integral equation, 
in general, are more efficient than the domain methods. As noted in [12], for N layers, the 
expression for the Green’s function would consist of a nested AM infinite series; hence, 
the evaluation of this Green’s function is somewhat computationally burdensome. 
Alternatively, the free-space Green’s function is used in [12] and [13] to avoid infinite 
series, but additional unknown charges on the dielectric interface and ground planes on top 
of the unknown charges on the conductor surface must be included, resulting in a larger 
moment matrix. Yet another approach to avoid an infinite summation is to solve the 
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integral equation in the spectral domain (SDA) [27], where the Green’s function is in a 
closed form; however, this approach can not be applied to conductors with finite 
thicknesses, such as MoL. 

In this thesis, the spatial closed-form Green’s function is used to avoid the 
evaluation of an infinite series without any additional unknowns used in the method based 
on the free-space Green’s function. A closed-form Green’s function for a multilayered 
dielectric medium was first introduced in [28]. This Green’s function utilizes a finite 
number of weighted complex images instead of an infinite number of real images required 
for the exact Green’s function. A closed-form Green’s function based on a finite number 
of weighted real images is first proposed in this thesis to avoid the use of expensive 
complex operations. 

1.4 Structure of the Thesis 

As mentioned in the previous section, the technique to solve an electrostatic 
problem (Laplace's equation) plays an important role in the quasi-TEM analysis, and 
methods based on an integral equation are used throughout this thesis to solve various 
electrostatic problems related to interconnections and discontinuities. All the integral 
equations encountered in this thesis are solved using the method of moments (MoM) [29], 
[30] with pulse basis functions and point matching (delta testing), and the moment matrices 
associated with the integral equations are constructed using an analytical formula for most 
cases, avoiding numerical integration or infinite summations. 

The core of an integral equation approach is the determination of the Green’s 
function. The exact Green’s function for a multilayered medium is often obtained by using 
the image theory, and it consists of an infinite number of images. Chapter 2 discusses an 
efficient expression of this Green’s function based on numerical approximation. This new 
expression of the approximate Green’s function uses only a finite number of images; 
hence, it is in a closed form. To obtain this closed-form expression, the spectral-domain 
Green’s function is first derived in this chapter; then, the spectral-domain Green’s function 
is approximated with real-valued exponential functions using the method based on the 
relaxation of curve fitting. The closed-form Green’s functions for a point, line, and semi- 
infinite line charges are then obtained by analytically converting the approximate spectral- 
domain Green’s function to the space domain. 
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Chapter 3 discusses the computation of the four parameters of a multiconductor 
transmission line, viz., the capacitance matrix C, the inductance matrix L, the resistance 
matrix R, and the conductance matrix G. The capacitance matrix C is computed from the 
free charge distribution on the surface of conductors, which is determined from the 
electrostatic analysis. To compute the inductance matrix L, an analogy between 

electrostatic and magnetostatic problems for uniform transmission line configurations is 
used. Hence, the conduction surface current distribution is computed by solving the 
equivalent electrostatic problem and then used to compute L. 

The resistance matrix R is also computed from the current distribution used in the 
computation of the inductance matrix by performing the perturbation analysis on this 
current distribution. Conventionally, the resistance matrix is defined in terms of power 
loss on conductors. The resulting matrix is nondiagonal in nature and is strongly 
dependent on the current excitations used in the computation. Thus, if the resistance matrix 
is obtained before the actual current distribution on the conductor has been determined, the 
result would not be too meaningful. In this chapter, the diagonal resistance matrix is 
defined in a manner such that it is relatively insensitive to the choice of current excitations 
compared to the nondiagonal resistance matrix, which is often computed using the 
perturbation on attenuation constants [31]. In addition to losses on the conductor traces, 
those due to imperfectly conducting ground planes are also incorporated into the resistance 
matrix. 


The remaining transmission parameter, the conductance matrix G, models dielectric 
losses and can be computed from the shunt current density. Since this current density is 
related to the normal component of the electric field at the surface of a conductor, which, in 
turn, is related to the surface charge density, the shunt current density can be obtained from 
the surface charge density of the lossless system when losses due to the imperfect dielectric 
media are small. 

Chapter 4 is devoted to modeling of various strip discontinuities. In particular, a 
method to compute the equivalent (excess) capacitance of junction discontinuities, such as 
open ends, step junctions, bends, and T-junctions, are presented. Unlike other 
approaches, only one integral equation is employed to handle the above discontinuities 
instead of formulating a different integral equation for each discontinuity type. The integral 
equation is formulated in terms of the excess charge distribution to avoid the numerical 
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instabilities associated with the total charge formulation, where the integral equation is 
formulated in terms of the total charge distribution [21]. 

Chapters 5 and 6 discuss the modeling of yet other discontinuity types, vias and 
crossovers, respectively. These discontinuities differ from the ones discussed in Chapter 4 
since conductor traces in these discontinuities could be located in the different dielectric 
layers. Furthermore, for a crossover case, traces are no longer electrically connected; 
hence, the equivalent capacitance of a crossover contains a mutual term in addition to two 
self-terms, and the coupled integral equations have to be solved instead of a single integral 
equation. Again, all integral equations are formulated in terms of the excess charge 
distributions. The utilization of the Fast Multipole Method (FMM) [32]-[36] in accelerating 
the MoM-based computation of the excess capacitance of a crossover is also considered in 
Chapter 6. 

Finally, the conclusions and some future work evolving from this thesis are 
presented in Chapter 7. 
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CHAPTER 2 

DERIVATION OF THE CLOSED-FORM GREEN’S FUNCTION FOR A 
MULTILAYERED DIELECTRIC MEDIUM 

2.1 Introduction 

Among the various electrostatic solution techniques mentioned in Section 1.3, an 
integral equation approach requires the most analytical effort, mainly due to the 
determination of the expression of the Green’s function for a multilayered dielectric 
medium. Fortunately, the spectral-domain expression of the Green’s function has already 
been found by several authors: the expression for full-wave analysis can be found in texts 
[1], [2], whereas the expression for a electrostatic problem can be found in several journal 
papers [3] -[6]. 

The conventional approach to obtain the expression of the Green’s function is the 
use of the Fourier transformation, in which the equation governing the physics of 
problems, the Helmholtz wave equation for full-wave analysis and the Laplace equation for 
electrostatic analysis, is converted to the spectral domain by transforming all but one of 
space variables. Then, the resulting equation, which is an ordinary differential equation in 
terms of the remaining one space variable, is analytically solved to obtain the expression of 
the Green’s function in the spectral domain. The major bottleneck of this approach is that 
the direct analytical inversion of this spectral-domain Green’s function to the space domain 
is often impossible. 

A simple but brute force approach for this inversion is the use of numerical 
integration [6], [7]. Although this approach is commonly used in full-wave analysis 
because of the complexity of the expression of the spectral-domain Green’s function [7], it 
is seldom used in static analysis since this approach is computationally intensive and does 
not allow an analytical expression for the Green’s function in the space domain. 1 Yet 
another simple approach is to expand the spectral-domain expression using the geometric 
series. Then, an analytical expression in the space domain is found by applying the inverse 
Fourier transformation formulas to the resulting series: the Sommerfeld identity (for 2-D 


'The comparison of the CPU time used in this approach and other alternative approaches is given 
in Section 3.2.2. 
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problems) and the Weyl identity (for 3-D problems) for full-wave analysis [Chapter 2, 2] 
and the equivalent identities ((2.16a) and (2.16b)) for static analysis. The resulting 
expression in the space domain can be shown to be identical to one obtained by applying 
the image theory directly in the space domain [8]. The major disadvantage of this approach 
or equivalently the approach based on the image theory is that the resulting expression 
involves the summation of a nested infinite series; for N layers, the expression would 
consist of a nested N - 1 infinite series as mentioned in Chapter 1. 

In this chapter, the closed-form expression of the Green s function in the space 
domain, which does not involve any numerical integration or nested infinite summations, is 
presented. The expression of the spectral-domain Green s function is first derived in the 
following section. This expression is different than the ones given in [3]-[6], and, as will 
be shown in Section 2.3, this form of the expression is more convenient for the purpose of 
obtaining the closed-form Green s function. Then, the spectral-domain Green s function is 
approximated with real-valued exponential functions, and the resulting approximate 
Green s function is analytically inverted to the space domain to obtain the closed-form 
Green s function. 

2.2 Derivation of the Spectral-Domain Greens function 

The cross-sectional view of the general geometry of a multilayered medium is 
shown in Fig. 2.1. An arbitrary number N d of nonmagnetic lossless dielectric layers are 
backed by two optional ground planes with possible top or bottom locations. All dielectric 
layers and ground planes are assumed to be infinite and uniform in the xz plane. 

Consider a unit point charge located at the mth layer at {x 0 , y 0 , z Q ) (Fig. 2.2). The 

nr\ 

three-dimensional Green s function G is the potential due to this point charge and 
satisfies the following Poisson s equation: 

V 2 G 3D (x, y, z\x 0 , y 0 , z 0 ) = ^y£(x - x 0 )S(y - y 0 )S(z - z Q ) (2. 1 ) 

on 

To assure the unique solution to the above equation, G has to satisfy the appropriate 
conditions at the boundary: G 3D is constant at the surface of the ground planes, and G 3D 
and the normal components of the displacement field must be continuous across the 
dielectric interfaces. Noting that the dielectric medium is uniform in two directions, we can 
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Optional Top Ground Plane 

y=d Nd 

e Nd* Mo Hi 

y= d Nd-i 



Optional Bottom Ground Plane 


Figure 2. 1 . The cross-sectional view of a multilayered dielectric medium. 



Figure 2.2. The geometric configuration used to determine the spectral-domain expression 
of the Green's function. 
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represent the Green's function and the point charge in the spectral domain in terms of their 
transforms in the x- and z-directions. The space-domain and the spectral-domain Green's 
functions are then related by 


+00+00 


G 3D (x,y f z\x 0 ,y 0 ,z 0 ) = -^2 J jdodpe- J ^ x - x ^ mz -^ ) G 3D (a,y,P^ 0 .y 0 ,z 0 ) (2.2a) 


+00+00 


G 3D (a,y,p\x 0 ,y 0 ,z 0 )= J jdxdze ja(x ■* o)+; ^ (z l ° ) G 3D (.x,y,dp 0 ,y 0 ,z 0 ) (2.2b) 


where G 3D (a, y, P\x a , y Q , z 0 ) is the 3-D spectral-domain Green's function and a and (3 are 
the transform variables associated with the x- and z-directions, respectively. Then, the 
corresponding equation for (2.1) in the spectral domain is written as 


(£_ 

K dy 2 


a 2 - (i 2 


G 3D (a,y,p\x 0 ,y 0 ,z 0 ) 


1 

£(y) 


S(y-y c ) 


(2.3) 


The general solution of the above equation is given by 


G 3D (Y,y\r 0 ) 


Ae ^ + Be ^ 

2e m y 



(2.4) 


where the first subscript m denotes the layer where the source is located, whereas the 
second subscript n will be used to denote the layer where the Green's function is evaluated. 
A and B are unknown expressions to be determined. Note that e m appears in (2.4) unlike 
(2.3), where e n appears. 

An identical expression can be obtained for the 2-D spectral-domain Green’s 
function G 2D (y,y\p 0 )by Fourier transforming G 2Z> (x, y|jc 0 ,_y 0 ) in the x-direction with the 
transform variable y assuming the layers are uniform in the jc-direction. Furthermore, since 
the unknown coefficients A and B are to be determined using the boundary conditions only 
in the y-direction, it is easily seen that the expressions for G and G must be identical 
under these Fourier transformations. Thus, in what follows, the superscripts for the 
spectral-domain Green's functions are omitted. 


i 
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Applying the boundary conditions at the dielectric interfaces and ground planes, 
(2.4) can be written as 


6(y, *,) = + r n ,^e-r^-y() 

IV 

(2.5a) 

6(7 , fc) = 

p? 

VI 

PN 

(2.5b) 

where 



n— 1 m 

A m ,n = \_ S j,j+i A m ,n = An,m J”J ^jJ-\ 

j=m j=n + 1 


(2.6) 

v+ _ T JJ+ 1 


(2.7a) 

r- _ T jJ~ 1 


(2.7b) 

f r JJ+t+ f j+yJ+2 e 2l,(d l- d ’+' ) 

w+, -. + r w+1 f J+u+2 ^-^> 


(2.8b) 

r w _, + r MJ _ 1 e 2nd '- 2 ~ dj -' ) 


(2.8b) 

e, - f/ _ 2e, 

hj £, + £y l,J £j + £j 


(2.9) 

Here, 7} j and T t j are the reflection and transmission coefficients. F,- , +1 takes the value 
of 0 or -1 if the yth layer is a half space, or the (/+l)th layer is a ground plane, respectively. 
f nn+ 1 is the generalized reflection coefficient, which is the ratio of the amplitudes of the 
potentials at y = d n due to the image charges located above and below y = d n . Al m and 


A m m are unknown expressions to be determined. The superscripts + and - are used to 
denote the cases for y>y 0 and y<y Q , respectively. 
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Using the facts, at y = y a , that ( 2 . 5 a) and ( 2 . 5 b) must be equal and that the normal 


component of the displacement field must be discontinuous by the magnitude of the charge, 
we can obtain expressions for A+ m and A~ m as follows: 

<m = M m [e +iy ° + f m ,m-\e +Y(2dm - 1 yo> ] 

(2.10a) 

Am,m = + f m>m+ 

(2.10b) 

where 


M =fl— r r 1 

m m 1 m,m- 1 J m,m+\ € j 

(2.11) 


The complete expression of the spectral-domain Green s function has now been 
derived. In ( 2 . 5 a), A+ „ can be physically interpreted as image charges (and the actual 
charge when m is equal to n) located below the observation point y, whereas the product of 
and f n n+ j can be interpreted as image charges located above the observation point. 
A similar interpretation can be given to ( 2 . 5 b). 


It is interesting to note that all of the above equations have the following form: 


fk(Y) = 


1 - c 2/*-l(7)c )t3 


( 2 . 12 ) 


where cj, ci, and C3 are some constants which satisfy 0 < C 2 fk-\(Y) e ^ 7 ' - 1. and the 
expression of /^_j(7) is in the same form as /*.(/) for k > 1 and is a constant function for 
k = 1 . The value of k depends on the number of dielectric layers; for instance, for N 
dielectric layers k takes values from 1 to N. Since 0 < C2fi c -\(Y) e ^ :j - 1 for all k, the 
geometric series can be used to expand /*(y), and the resulting series is a nested N - 1 
infinite series. The entire expression of the spectral-domain Green s function can then be 
written in terms of this series, and each term of the series can be conveniently inverted to 
the space domain using the integration formulas ( 2 . 16 a) and ( 2 . 16 b) to be given in the 
following section. It can further be shown that the resulting expression of the Green s 
function in the space domain is identical to the one obtained by using the image theory. 
Hence, the exact expression of the Green s function in the space domain generally consists 
of a nested infinite series. 
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In the following section, the closed-form spatial Green s function, which avoids 
this nested infinite series, is derived based on the expression of the spectral-domain 
Green s function obtained in this section. 

2 . 3 Derivation of the Space-Domain Green s function 

In this section, the closed-form expression of the Green s function in the space 
domain is obtained by applying the exponential approximation to the spectral-domain 
Green s function. To obtain the closed-form expression of the Green s function in the 
space domain, the spectral-domain expression of the Green s function derived in the 
previous section is first rearranged by factoring out all y and y a dependencies as follows: 

G(y, y\r Q ) = — — ( Ki (y, m, n)e Y(y+y ° ~ 2 ‘ dn ) + K \ (y, m, n)e Yiy ~ y ° +2(dm -> ~ d " )} 

2 £ m y x 

■¥K^{y,m,n)e Y( '~ y+y °^ + K^{y,m,n)e Y( '~ y ~ y ° +2dm ~ x ^ y>y 0 (2.13a) 

G(y, y\r Q ) = — — (iff (y, m, n)e Y(y+y ° ~ 2dm) + h :j(y, m, n)e Yiy ~ y ° ) 

2 e m y ' 

+K$ (y, m, n)e Y( '~ y+yo+2( ' dn ~ 1 ~ dm ^ + K^ t n 4 e Y( '~ y ~ y ° +2dn ~^'j y < y G (2.13a) 

where 

n - 1 

^4 ( Y> m > n ) ~ 

j=m 

n- 1 

(7> m > n ) = ^m^n,n+l^m t m-l J J^jj+1 

j=m 

n— 1 

K${y,m,n) = M m Y\Sj, j+ \ 

j-m 

n-l 

K 4 (Y» m > n ) = f m m-1 | Sj, j + 1 

j-m 


(2.14a) 
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m 

K\ (7, m t ft) = M m f'm f m+\ ik-. 

j=n+\ 

m 

K2(y,m,n) = M m 

;=n+ 1 

m 

K3 (/> m> n) = M m f mm +\f nn _\ n«H 

j=n + 1 

m 

KHy,m,n) = M m f„^Y\. S i.M < 214b > 

y=«+i 


The determination of the closed-form spatial Green's function can now be preceded by 
approximating the above four coefficient functions K.f(y,m,n) using exponential 
functions. It is important to mention that although K.f{y,m,n ) is dependent on m and n, it 
is not a function of the source and observation locations, y and y 0 ; hence, the 
approximation can be performed without any prior knowledge of the geometry of the 
conductors. 

One physically intuitive approach to approximate the potential due to a charge in the 
layered medium may be the use of a finite number of weighted image charges in the 
homogenous medium, which is equivalent to approximating the coefficient functions 
K~{y, m,n) with exponential functions. These weighted images can be either complex or 
real depending on whether complex-valued or real-valued exponential functions are used in 
the approximation. The equivalence between the weighted image charges in the space 
domain and exponential functions in the spectral domain will be shown later in (2. 18a) and 
(2.18b). 

In electromagnetic analyses, the complex-valued exponential functions are often 
used for pole-zero modeling of signals, such as an electromagnetic-scatterer response. The 
least-square formulation of this exponential approximation results in nonlinear equations 
and can only be solved by iterative methods, such as gradient descent procedures or the 
Newton method. Due to the computational inefficiency of these algorithms, some other 
suboptimal noniterative techniques are proposed: the least-squares Prony method and the 
generalized pencil-of-fiinction (GPOF) method [9], [10]. These suboptimal methods are 
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used to obtain the closed-form Greens function for full-wave analysis [11]-[13] 2 and 
further applied to obtain the closed-form Greens function for electrostatic analysis [8]. 
Although these algorithms are noniterative, their computation involves matrix inversions 
and a polynomial factoring or a solution of the generalized eigenvalue problem, which still 
can be considered as computationally inefficient. 

Fortunately, the four coefficient functions Kf{y,m,ri) in (2.13a) and (2.13b) are 
nonoscillatory and smooth functions of y; hence, each coefficient function can be 
sufficiently approximated with real-valued exponential functions instead of complex- valued 
exponential functions, avoiding computationally expensive complex operations. The real- 
valued exponential approximation method described in [14] is employed in this thesis. The 
method is based on the relaxation of curve fitting, and the details of the procedure are given 
in Appendix A. Although this method is simple and iterative in nature, it converges to 
reasonable accuracy in a few iterations and requires much less computation time as 
compared to those for the previously mentioned methods. 

It can be seen that a pole exists at y = 0 for a medium with both top and bottom 
ground planes, and K.f{y, m, n ) can no longer be accurately approximated with exponential 
functions. Thus, a special treatment is required for this case to extract the pole from 
Kf{y, m, n). In the following subsection, the closed-form Green s function is obtained for 
cases with no ground planes or only the bottom ground plane, and the subsequent 
subsection discusses the derivation of the Green s function for a case with both top and 
bottom ground planes. 

2.3.1 Closed-form Greens functions for geometries without any ground 
planes or with only the bottom ground plane 

When the dielectric layers are not backed by both top and bottom ground planes, the 
four coefficient functions K.f(y,m,n) do not have any poles; furthermore, they are 
nonoscillatory and smooth functions of y . Hence, they can be approximated using the 
exponential approximation method discussed in Appendix A as follows: 


^Strictly speaking, the closed-form Green s function does not exist for a full-wave case since y and 
y 0 dependencies cannot be removed from the coefficient functions before the approximation. 
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Kf(m, n,y) 



f^'j ,, a m,njY 

m.n.r 


(2.15) 


where N± ni denotes the number of exponential functions used in the approximation of 
K±(y, m,n), which typically ranges from 5 to 10. In general, K±(y, m,n) have asymptotic 
values, and an analytical extraction of these values should be performed to increase the 
accuracy of the approximation or to reduce the computation time. These asymptotic values 
can be easily obtained either analytically or numerically. Furthermore, some of these 
coefficient functions are often zero or one, and these properties can also be explored for 
further computer time savings. 


Once the exponential approximation is performed in the spectral domain, the closed- 
form Green s functions for 2-D and 3-D can be obtained in the space domain by using the 
following inverse Fourier transformation formulas: 


+°o 


-ln(p) = -ln(V*W) = ^ jrfje-'’* 

— oo 



1 _ 1 

"r = ^ x 2 +y 2 +z 2 


+00 +°° — ylyl 

^ J J dadPe- j(ocx+lk) ~- 

— oo — oo 


(2.16a) 


(2.16b) 


where y = ^a 2 + ft 2 in (2.16b). The above identities can be easily derived by 
considering the potentials due to the unit point and line charges in a homogeneous medium 
without any ground planes, and the corresponding identities for a full-wave case are 
Sommerfeld s and Weyl s identities, respectively. Applying the above formulas to (2.13a) 
and (2.13b), the closed-form 2-D and 3-D Green s functions are written as 


G 2 D (P \Po) = - 




(2.17a) 


(2.17b) 


For i=l and y > y 0 , the expressions of 


f 2D,± (r\r 0 ) and f? D,± (r\r 0 ) are given by 
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r2D,+ 

J 1 




(p\po ) = X 1 • (x -x 0 ) 2 +iy + y 0 - 2d n + a^ n j Y 

7=1 


(2.18a) 


j m ’ W ’ 1 + , 2 (2.18b) 

,=1 ^(x-x 0 ) z +(y + y 0 - 2 d n + a^ nl ) 1 + (z - z 0 ) 


Similar expressions can be obtained for 


ff D,± (r\r 0 ) and f? D,± (r\r 0 ) for other values of i. 


The derivations of the closed-form Green s functions for a point charge and a line 
charge are now completed for geometries without both top and bottom ground planes. 
Considering the forms of (2.18a) and (2.18b) it is clear that the exponential functions used 
to approximate the Green s function in the spectral domain correspond to the weighted 
images in the space domain. 


In the computation of the equivalent capacitances of interconnection discontinuities, 
. the Green s function for a uniform semi-infinite line charge, G semt (r\r 0 ,^) , is required to 
formulate the integral equation in terms of the excess charge distribution. The closed-form 
Green s function for a semi-infinite line charge is derived in the rest of this subsection. 

To derive G semi (r]r ot %), the auxiliary Greens function for a line charge with 
polarity reversal is employed [15]. Consider a uniform line charge, which starts from 
z = € and is infinitely extended in the positive z-direction (see Fig. 2.3); then 
G semi ( r \r 0 ,%) can be expressed as 

G” mi fr|r D .|;= i[G 2c fp| p„)+ C '7r|r„,^] (2.19) 


where G p (r\r 0 ,£) is the Green s function for a line charge with an abrupt polarity reversed 
from minus to plus at z = £. Since the closed-form expression of G 1D (p\p Q ) can be 
obtained using the previous technique, the closed-form expression of G semt (r\r 0 ,%) can be 
determined once the closed-form expression is derived for G p {r\r 0 ,E l ). The expression for 
G p (r\r 0> 4) is rather easily obtained by integrating the potential due to a point charge: 
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(a) 


(b) 


Figure 2.3. The decomposition of (a) the uniform semi-infinite line charge to (b) an infinite 
uniform line charge and an infinite line charge with the polarity reversal at 

z = l- 

^ 4 

G<’(r\r 0 ^) = - J C 30 (r|r 0 ) + Jg 3d (^) = • ± (Hr<,,|) (2-20) 

-oo £ m 1=1 

The integration can be performed using the following formula [(3.3.40), 16]: 


f , £ - = In 

x + ^Jx 2 ±a 2 

J V * 2 ±a 2 



( 2 . 21 ) 


Again, for i = 1 and y>y 0 , / ( p ’ + (r|r 0 ,^> is given by 


<», i 

(+•<,,£) = X c mir ln 

;=i 


^(jr-jt 0 ) 2 +(y + y 0 -2d„+a^ 1 ) 2 +(z-^) 2 + (z-j) 

(* - j: 0 ) 2 + (y + y 0 ~ 2d n + a mii,l ) 2 + (* - £) 2 - ~ <a) 

( 2 . 22 ) 


2.3.2 Closed-form Greens functions for geometries with both top and 
bottom ground planes 


As mentioned earlier, when both top and bottom ground planes are present, all of 
the four coefficient functions K+(y, m, n ) are still nonoscillatory but contain a pole at y= 0. 
Since exponential functions are bounded, they cannot be used to approximate unbounded 
functions; therefore, this singularity must be extracted prior to the exponential approxima- 
tion to preserve good approximation results. In [11], the complex exponential approxima- 
tion has been performed without considering this singularity and used to obtain the closed- 
form Green s function. As a consequence, the resulting capacitance values from this 
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Green s function were shown to contain large errors. The extraction of this pole can be 
accomplished by rewriting G(y,y\r 0 ) in the following manner: 

G(y, y|r 0 ) = R m>n G h (y, y\r 0 ) + G (y, y\r 0 ) (2.23) 

where G h (y,y\r 0 ) is the spectral-domain Greens function for a homogeneous medium 
with the same ground planes, i.e., all dielectric layers are replaced by the source layer. 
Again, G h (y,^r 0 ) contains a pole at y= 0. R m n is a constant which is determined such 
that G (y,y|r G ) is a well-behaved function without any poles. R m n can be obtained either 
numerically or analytically by taking limits of G(y,y\r 0 ) and G h (y, y\r 0 ) as y — » 0. 


Now the technique used in the previous section can be applied to obtain the closed- 
form expression for G (y,y|r c ) in the space domain, and the expressions of G(y,y\r 0 ) are 
obtained once the space-domain expressions of G h (y,y\r a ) are determined. Since the 
medium is homogeneous for G h {y,y\r 0 ), the expressions can be easily obtained using the 
image theory and are given by 


0 2D \ P \Po)-~ t 


In 




k=- 


t](x-x 0 ) 2 +(y-y 0 -2kh) 2 
^/(x-x 0 ) 2 +(y + y 0 -2 kh) 2 


(2.24a) 


^(x - x 0 ) 2 + (y - y a - 2khf + (z - z 0 ) 2 


^/(x-x 0 ) 2 +(y + y D -2k/i) 2 -l-(z-z f> ) 2 


(2.24b) 


k=~<> 


In 


j(x-x 0 ) 2 +(y-y o - 2kh) 2 +(z- £) 2 +(z-$) 
J(x - x 0 ) 2 + (y - y 0 - 2Mi) 2 + (z - £) 2 -(z-£) 


^(x-x 0 ) 2 +(y + y 0 -2kh) 2 +(z-%f -(z-Z) 
A /(x-x 0 ) 2 +(y + y 0 -2k/i) 2 +(z-^) 2 + (z -£) } 


(2.24c) 


Unfortunately, all expressions are written in terms of infinite series. G 2D ' h {p\p Q ) can be 
alternatively expressed using a closed-form formula [17], but this closed-form expression 
requires numerical integration when the moment matrix is computed, unlike (2.24a), which 
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can be integrated analytically (see Section 2.4). The detailed discussion of alternative 
expressions for G 2D,h (p\p 0 ) is given in Appendix B. For the remaining Green s 
functions, G 3D,h (r\r 0 ) and G p ' h (r\r 0 ,^), such a closed-form formula does not exist at least 
without special functions and an infinite-series expression cannot be avoided when both top 
and bottom ground planes are present. However, the expressions for G 2D (p\p 0 ), 
G 3D (r \r D ), and G p (r\r 0 ,%) given in this subsection are still numerically more efficient than 
the ones obtained from the conventional image method since a nested infinite series ( N - 1 
nested infinite series for N layers) of the conventional method is reduced to a simple infinite 
series without any nesting as shown in the above equations. For this reason we shall still 
refer to the Green s functions expressed by (2.23), (2.24a), (2.24b) and (2.24c) as closed- 
form Green s functions. 

2.4 Closed-form Integration Formulas for the Elements of Moment 

Matrices 

When integral equations are solved using the method of moments (MoM), the 
elements of the moment matrix are computed by integrating the Green s function over basis 
and testing functions; the closed-form formula for this integration is discussed in this 
section. The general forms of the integrations required to construct the moment matrix are 

J J T(p)G 2D (p\p s )B(p s )dp s dp (2.25a) 

l,h 

f \T(r)G 3D (r\r s )B(r s )dr s dr (2.25b) 

44 

where T and B are testing and basis functions, l s is the source line segment, and l t is the 
testing line segment where the potential is evaluated for 2-D problems. Similarly, Q s and 
Q t are the source and testing patches for 3-D problems. The method of collocation uses 
the delta testing function (point matching) to reduce the above double integrations to a 
single integration; it is used in this thesis whenever the moment method is employed. The 
testing point is usually chosen at the center of a segment or a patch. 

For basis functions used to expand the surface charge density, pulse-type 
functions, which are commonly used in the method of collocation, are used in this thesis. 
Then, (2.25 a) and (2.25b) simplify to 
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jG 2D (p c \p s )dp s (2.26a) 

h 

J G 2D {r c \r s )dr s (2.26b) 

O t 


where p c and r c are the centers of the testing segment and patch. Similarly, the integral 
associated with G semi (r c \r s ,%) can be written as 


fG semi (r c \r s ,4)dp s (2.26c) 

h 


Equation (2.26a) appears in the computation of the moment matrix for 2-D problems or the 
right-hand side vector (excitation vector) computation for the equivalent capacitance of a 
crossover. Equation (2.26b) only appears in the computation of the moment matrix for 3-D 
problems, whereas (2.26c) only appears in the computation of the right-hand side vector 
for the equivalent capacitances for various junction discontinuities. 


After substituting (2. 17a) and (2. 18a) into (2.26a), the line integrals associated with 
each term in the summation can be put into the following form with manipulations in y, y 0 , 
and the terms due to the exponential approximation: 

J = J ln(|p - p’|)<tf = J ln(P)dr (2.27a) 

c c c 


Similarly, for 3-D problems we have 


JJ 


d£ 

^l(x-x ) 2 +(y-y') 2 +(z-z') 2 



(2.27b) 


The evaluations of (2.19a) and (2.19b) over an arbitrary polygonal patch and a line 
segment are well-known and the closed-form formulas are given in [18]. On the other 
hand, (2.26c) can also be integrated analytically using the following integration formula 
[19]: 
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‘2 

Jl" 


Va 2 +fc 2 +/ 2 +a 
Va 2 +6 2 +/ 2 -a 


k = 2W 


tan 


-l 


al i 


a 2 4- f> 2 + /j 2 


-tan 


-1 


a/o 


mV* 


2 +b 2 +l 2 2 j 


( 


+ 2aln 


h + i 

la 2 + b 2 + 1 2 2 


la 2 + b 2 + 1 2 


+h^ n 



a 2 +b 2 +l 2 2 + a 

-l jin 


a 2 + b 2 + l x 2 + a ^ 

vAl 

l a 2 +b 2 +l 2 2 -aj 

J 

a 2 +b 2 + 1 2 -dj 


(2.27c) 


Finally, using the above integration formulas all elements of the moment matrices 
and the right-hand side vectors encountered in this thesis can be performed analytically. 

2.5 The Comparison of the Exact and Approximate Greens Functions 

In this section, the comparison of the closed-form Green s function, which is 
approximate, and the exact Green s function is presented. As a first example, the medium 
shown in Fig. 2.4 is considered. The maximum number of 1 1 exponential functions are 



Ground Plane 


Figure 2.4. The first geometry used to test the closed-form Green s function. 
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used in the approximation. The observation and source points, y and y 0 , were taken to be 
0.6 mm and 1.6 mm, respectively. Figure 2.5(a) shows the comparison of the 
approximate and the exact Green s functions, and Fig. 2.5(b) shows the relative errors; the 
maximum relative error was less than 0.023 %. It is important to observe that although the 
exponential approximation might fail for the large argument case due to its fast decaying 
nature, by extracting the asymptotic values and the y and y 0 related exponential factors 
from the coefficient functions, the limiting behavior of the overall approximated Green's 
function would still remain accurate for the large values of y. 


As a second example, the two-layer medium case shown in Fig. 2.6 is considered. 
A maximum number of exponential functions used in the approximation was five, y and y 0 
were again taken to be 0.6 mm and 1.6 mm. The comparison is shown in Figs. 2.7(a) and 
2.7(b), and an excellent agreement was found. The approximate closed-form Green s 
function is also compared with the exact Green s function in the space domain. The exact 
2-D Green s function is obtained using the image method and is given by 


G 2D (p\p 0 ) 


(1 + n y J {x-x o ) 2 +{y-y 0 -(2k-l)h) 2 ' 
4n£ o £ {(x-x 0 ) 2 +(y-y 0 -(2k-3)hj 1 j 


y>h 


(2.28a) 


G 2D (p\p 0 ) 


(i+n 
4 ne a 



(*-*o ) 2 

(*-*o ) 2 


+ (y + (2fc-l)/i) 2 > 
+ (y + (2k-3)h) 2 j 


y<h 


( 2 . 28 b) 


where 


r = £2 — £L ( 2 . 28 c) 

e o + e l 

In the above expressions, it is assumed that the line source is located at the dielectric 
interface. 

The comparison results are shown in Figs. 2.8(a) and 2.8(b), and the results were 
in good agreement. This example justifies the validity of the approximation of the Green s 
function in the spectral domain. It is observed from numerous approximations that the 
smaller number of dielectric layers required fewer exponential functions to approximate, as 
expected. In fact, for one layer with a ground plane, the closed-form Green s function 
becomes exact since there is only one image charge. 


The Relative Error Potential (V) 
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Comparison of the Green's functions 



(a) 


The Relative Error in the Approximation 



(b) 


Figure 2.5. (a) Comparison of approximated and exact Green’s functions in the spectral 
domain for the medium shown in Fig. 2.4 and (b) the plot of the relative error. 
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Ground Plane 


Figure 2.6. The second geometry used to test the closed-form Green s function. 




The Relative Error Potential (V) 
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Comparison of the Green's Funtions 



(a) 

The Relative Error in the Approximation 



Figure 2.7. (a) Comparison of approximated and exact Green's functions in the spectral 
domain for the medium shown in Fig. 2.6 and (b) the plot of the relative error. 




The Relative Error 
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Comparison of the Green's Functions 



The Relative Error in the Approximation 



Figure 2.8. (a) Comparison of approximated and exact Green’s functions in the space 
domain for the medium shown in Fig. 2.6 and (b) the plot of the relative error. 
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For the final example, the stripline case shown in Fig. 2.9 is considered. In 
general, a slightly larger number of exponential functions was required for a strip case, 
where both the top and bottom ground planes were present, and a maximum number of 13 
exponential functions were used in this approximation. Both y and y 0 were taken to be 
0.1 mm. To compare the approximation results, G (y,yjr 0 )in (2.23) is considered instead 
of G(y,y\r 0 ). The comparison results are shown in Figs. 2.10(a) and 2.10(b), and a good 
agreement was found. 

For most numerical examples given in this thesis, the exponential approximation 
took less than a second of the CPU time in Sparc 10 and Alpha workstations; thus, the 
CPU time used in the approximation was negligible compared to the one for the 
construction or the inversion of the moment matrix. 

Finally, it should be mentioned that since the method of moments is an approximate 
way of solving an integral equation, only a moderate accuracy is required for the Green s 
function. In general, five exponential functions for cases with the presence of either one of 
the top and bottom ground planes and seven exponential functions for cases with the 
presence of both the top and bottom ground planes were enough to solve most capacitance 
problems. 


y t 


I Ground Plane 



Ground Plane 


Figure 2.9. The second geometry used to test the closed-form Green s function. 





Comparison of the Green's Functions in 
the Spectral Domain 



Figure 2.10. Comparison of approximated and exact Green's functions in the spectral 
domain for the medium shown in Fig. 2.9. 
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2.6 Physical Interpretation of the Closed-Form Greens Functions 

As shown in the previous section, the closed-form Green s function well 
approximated the exact Green s function, which requires an infinite number of images, 
using only a finite number of weighted images. In this section, the underlying physics of 
this closed-form Green s function is discussed to understand the difference between the 
images due to the image principle and the exponential approximation. 

Let us first assume that we are only interested in the circular region around a point 
source as depicted in Fig. 2.11(a). Then, the infinite number of images in the exact 
Green s function may be truncated to a finite number of images to evaluate the potential due 
to this point source inside of the region of interest up to a certain desired accuracy as shown 
in Fig. 2.11(b). This finite number of images is spaced and weighted such that its 
extension to the infinite number of images satisfies the boundary condition at the dielectric 
interface over infinite distances. The boundary condition on the portion of the dielectric 
interface that lies inside the region of interest may still be satisfied by employing a new set 
of images, which are nonuniformly spaced and weighted. Because of the nonuniform 
spacing and weighting, the smaller number of images than the one for Fig. 2. 11(b) may be 
required to satisfy the boundary condition as illustrated in Fig. 2.11(c). Thus, the finite 
weighted images in the closed-form Green s function may be used to evaluate the potential 
distribution inside this region of interest. 

Since the closed-form Green s function utilizes only a finite number of images, 
these images are located over a finite region; thus, these images cannot satisfy the boundary 
condition at large distances and may not be used to evaluate the potential at such large 
distances. The similar situation applies to the truncated version of an infinite number of 
images. However, the valid range of the weighted images resulting from the spectral- 
domain exponential approximation turns out to be large enough for most practical 
problems. For instance, the mutual capacitance between the two conductors computed 
using the closed-form Green s function is accurate until the separation distance between the 
two conductors becomes large enough that the mutual capacitance becomes negligible as 
compared to the self-capacitances. 
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(c) 

Figure 2.1 1. (a) A point source in inhomogenous media and its equivalent systems using 
(b) the image theory and (c) weighted images. 
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2.7 Summary 

Numerically efficient forms of the Green’s functions for a point charge, a line 
charge, and a semi-infinite line charge embedded in a multilayered dielectric media with two 
optional ground planes were presented in this chapter. The presented Green’s functions are 
approximate and utilize a finite number of weighted images instead of an infinite number of 
images used in the conventional exact Green’s function. The analytical integration 
formulas to integrate these Green’s functions to evaluate the elements of the moment 
matrices and the excitation vectors were also presented. 
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CHAPTER 3 

COMPUTATION OF THE TRANSMISSION PARAMETERS OF A 
MULTICONDUCTOR SYSTEM 

3.1 Introduction 

As the word “transmission” implies, electromagnetic fields associated with a 
transmission line are dynamic in nature. However, for TEM lines, the transverse 
distribution of the fields at any instant of time is identical to that for the static solution. As a 
consequence, the four parameters for multiconductor TEM transmission lines, viz., the 
resistance matrix R, the inductance matrix L, the conductance matrix G, and the 
capacitance matrix C, may be derived from a static analysis with good accuracy. Similarly, 
under the quasi-TEM approximation, the spatial distribution of the fields in a multilayered 
dielectric medium is essentially identical to that predicted by the static analysis; hence, the 
R, L, G, and C matrices obtained from the static analysis still represent good 
approximations to the quasi-TEM transmission line parameters. 

It is well-known that by making an analogy between electrostatic and magnetostatic 
problems for uniform transmission line configurations, both the charge and current 
distributions can be obtained from the electrostatic analysis [1], and the capacitance matrix 
C and the inductance matrix L follow from these distributions. Among various methods to 
determine the charge distribution on a multiconductor system [l]-[9], the most common 
method is the integral equation approach. A brief comparison of various methods used for 
electrostatic problems was presented in Section 1.3. In this chapter, an integral equation is 
formulated in terms of the closed-form Green’s function derived in the previous chapter; 
thus, the charge distribution is determined in an efficient manner such that no nested infinite 
summations nor numerical integrations are performed. To determine the current 
distribution, the equivalent electrostatic problem is solved using the same technique 
employed in the computation of the charge distribution. 

Losses due to imperfect dielectrics by themselves do not alter the quasi-TEM nature 
of multiconductor transmission lines; hence, the conductance matrix G can still be 
computed from the solution of the electrostatic problem by introducing a complex dielectric 
constant to account for the finite loss tangent [10], [11]. This procedure, however, is 
computationally inefficient (compared to the proposed method in this chapter) as it involves 
complex operations. An alternate approach, which is based on the conductance analogy of 
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capacitances, is often applied to single lines (two conductors) embedded in a homogeneous 
medium. In this chapter, the above concept is generalized to compute the conductance 
matrix G of n-lines embedded in a multilayered dielectric medium. 

Unlike the case for the computation of the conductance matrix, where the 
electrostatic problem could have been formulated with a complex permittivity, a formulation 
for the resistance matrix R in the quasi-static regime is not so evident since the quasi-TEM 
approximation inevitably neglects the longitudinal components of fields, which must be 
present due to the conductor loss. The incremental induction method [12] is commonly 
used to compute the resistance for single line cases; however, it cannot be applied to 
general multiconductor cases. Yet another commonly used method is one based on the 
perturbational analysis on attenuation constants [10]. This approach first computes the 
attenuation constant by taking the ratio of the perturbed modal power loss and the total 
modal power of a lossless system; then, the resistance matrix is obtained by solving a set of 
N 2 linear equations for N modes, and the resulting resistance matrix is shown to be 
nondiagonal. Although this approach is variational (second-order accurate) in terms of 
attenuation constants, it is no longer variational for the resistance matrix. 

As indicated in [1 1] and Section 3.5 in this chapter, the resistance matrix of coupled 
transmission lines is a nondiagonal matrix and is strongly dependent upon the choice of the 
current excitations used in the computation; the value of the resistance matrix varies 
significantly as different current excitations are used. This undesirable phenomenon is 
more prominent for the off-diagonal elements of the resistance matrix, and prevents the 
computation of R without the knowledge of the actual current distribution, which is known 
only in the simulation time. Consequently, the nondiagonal matrix form of the R matrix 
may be of no practical use. In this chapter, the diagonal matrix form of the resistance 
matrix, which is computed from the total power loss, is proposed. The resulting matrix is 
shown to be relatively insensitive to the choice of current excitations; hence, the proposed 
diagonal matrix is more suitable for the computation of the resistance matrix. Furthermore, 
the procedure to compute this matrix does not involve eigenvalue analysis unlike the 
perturbational approach mentioned previously. 

In the following four sections, the details of the previously proposed methods to 
compute C, L, G, and R matrices of a multiconductor system is presented. Figure 3.1 
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Optional Top Ground Plane 


y =d Nd 



Optional Bottom Ground Plane 

Figure 3.1. The cross-sectional view of a possible configuration of multiconductors 
multilayered dielectric medium. 


in a 


demonstrates the general geometry of a multiconductor system embedded in a multilayered 
dielectric medium. The dielectric layers are the same as the ones considered in the previous 
chapter. An arbitrary number N c of conductors are placed throughout the layers, and the 
cross sections and planar geometries of the conductors can be arbitrary as long as their 
boundaries can be described with a piecewise linear function. 

3.2 Computation of the Capacitance Matrix 

3.2.1 Theory 

An impressed potential on conductors results in free charge accumulation on the 
surfaces of conductors, and the electrostatic potential 0(r) at any point except inside the 
conductors is then related to this surface charge density q(p) per unit length via the 
following integral equation: 

0(r)= Jg 2D (p|p' W) dr' = (g 2D ,<?) (3-1) 

a 

where Q denotes the contours of all conductors except for ground planes. G 2D (p\p ' ) is 
the 2-D closed-form Green’s function for a multilayered medium discussed in the previous 
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chapter, and it accounts for polarization charges on the dielectric and conductor interfaces 
and free charges on the surfaces of ground planes. The integration is symbolically written 
as (•,•) to simplify the notation. To solve the above integral equation numerically, the 
method of collocation is used. First, the conductors are approximated with polygons and 
the unknown charge density is expanded with the pulse basis functions; then, the point 
matching technique is applied to the centers of each basis functions. The integral equation 
(3.1) can now be put into the following matrix form: 


where 




’ q\ " 


M l,l M l,2 * M 1 ,N c 

’ ^1 ' 

v 2 

= M 

h 



M 2 ,i • 

h 

• 


• 


• 

• 



i 

£ 

i 


M Nc ,t M N c ,N c _ 

' ^1 

£ 

1 


(3.2a) 


Vi=[V h - . V,f (3.2b) 

Qi={<ll • • . (3.2c) 

= jG 2D (pi’ p |p)4t> (3. 2d) 

r i 


where the superscript T denotes the transpose, Vj and q t are vectors of size N it and Mj j 
is an Ni by Nj matrix, where N ( is the number of basis functions used to discretize the rth 
conductor. V i is the voltage of the rth conductor with respect to the ground planes, q l j is 
the unknown coefficient associated with the y'th basis function of the rth conductor, W is 
the q\h line segment of the y'th conductor, and p l c ' p is the center point of the pth basis 
function of the rth conductor. The expression of G 2D (p\p') is given in Section 2.3, and 
the closed-form formula for the integration in (3.2d) is discussed in Section 2.4. Now, 
given the excitation voltages, the corresponding charge distribution can be determined by 
solving the linear system of Equations (3.2). 


The capacitance matrix C is defined to relate the total free charges on the conductors 
to the voltages in the following manner: 


CV = Q 


(3.3a) 


where the ith elements of V and Q are the voltage and the total charge of the fth conductor, 
and the ith element of Q, Q, is given by 

Ni 

a = < 3 - 3b) 

Here, /j and q l j are the length and the charge density coefficient of the yth segment of the 
ith conductor. 

The capacitance matrix can now be determined by solving the charge distributions 
for N independent voltage excitation vectors and can be put into the following matrix form: 

C = QV -1 (3.4) 

where the ith columns of V and Q are the ith voltage excitation vector and charge vector. 
To avoid the matrix inversion, the identity matrix is chosen for the V matrix in this chapter. 
With this choice of V matrix, C is simply equal to Q, and the following physical properties 
of the capacitance matrix can be deduced: 

■ Since Q is symmetric due to reciprocity, C is also symmetric: 

C,j = Cjj (3.5a) 

■ Since the diagonal elements of Q are the charges on the excited 
conductors and the off-diagonal elements of Q are the induced charges 
on the rested conductors, the diagonal elements of Q are positive and all 
of the off-diagonal elements are negative: 

C, ,• > 0 and Qj < 0 for j * i (3.5b) 


Furthermore, the magnitude of the sum of the induced charges must be 
smaller than the charges on the excited conductor; hence, C is a positive 
definite matrix: 
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The capacitance matrix defined by (3.3) is called the Maxwellian capacitance, and 
the physical self- and mutual capacitances are related to this Maxwellian capacitance matrix 
by 


<=i.i — 

(3.6a) 


(3.6b) 

J*i 



So far, we have implicitly assumed that there is at least one ground plane, which is 
the reference conductor. If there are no ground planes, any one of the conductors can be 
chosen as a reference conductor, and the terms corresponding to this reference conductor 
should be removed in (3.3a) and (3.4). Unfortunately, when there are no ground planes, 
the physical condition in which the total charges in the system must add to zero for 2-D 
problems is not necessarily satisfied by the integral equation (3.1), which was naturally 
enforced by the Green’s function when there is at least one ground plane. Thus, this 
condition must be incorporated into (3.1) to obtain the correct charge distribution. To 
enforce this condition, another row is added to (3.2) as follows: 




• • • • 

• 

• 


• • • • 

• 

• 


• • • • 

• 

0 


1 1 1 0_ 

1 

1 

1 


where V re j is the potential value at the reference conductor, which is yet to be determined. 

It is interesting to note why the above physical condition has to be satisfied in 
addition to the boundary condition embedded in (3.2) in spite of the facts that all necessary 
boundary conditions appeared to be embedded in (3.2) and that the solution of the T ^p1a<r 
equation is unique once the boundary conditions are satisfied. The answer to this question 
is that since the magnitudes of the fields have to be finite everywhere except at the source 
region, they have to be finite at infinity also, and this boundary condition at infinity turns 
out to be missing in the integral equation (3.1) for cases in which no ground plane is 
presented. This boundary condition at infinity can be enforced by setting the total charge in 
the system to zero. It should be noted that when there is one ground plane, this boundary 
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is embedded in the Green’s function, and when dielectric layers are bounded by two 
ground planes, this boundary condition no longer has to be satisfied. 

3.2.2 Numerical examples 

In this section, various numerical examples arc given to verify the method 
discussed in the previous section. As a first example, a single microstrip line shown in 
Fig. 3.2 is considered. The line thickness was taken to be infinitely thin. The capacitance 
value is computed with various numbers of basis functions, and the result is compared with 
those for the iterative spectral-domain technique [13], [14] in Fig. 3.3. This spectral 
domain method solves the integral equation iteratively in conjunction with the minimization 
in the boundary condition error. The maximum number of exponential functions used in 
the approximation of each coefficient function was seven. The charge density on the 
microstrip is plotted in Fig. 3.4. 

A more complex geometry, a three-conductor system in a layered medium, shown 
in Fig. 3.5, is considered. The number of basis functions used was fifty for each 
conductor. Again, the maximum number of exponential functions used in the 
approximation was seven. A comparison with the results obtained from [4] is shown in 
Table 3.1. In [4], the spectral-domain Green’s function is numerically integrated to obtain 
the space-domain Green’s function using a Gaussian quadrature formula in conjunction 
with an analytical asymptotic extraction. The potential distribution in the system with the 
excitation at the center conductor is plotted in Fig. 3.6. The white spaces in Fig 3.6 
represent the conductors. 

A ten-conductor transmission-line system above a thick dielectric substrate, shown in Fig. 
3.7, is also considered. The total number of 300 basis functions was used to represent the 
unknown charges, and nine exponential functions were used to approximate each 
coefficient function of the Green's function. Table 3.2 shows the computed results. The 
same structure is considered in [3] using the free-space Green’s function with the basis 
functions, which incorporate the edge singularities of the charge near the comers of the 
conductors. In [3], data are obtained using 160 and 190 basis functions for conductors and 
dielectric interfaces, respectively. The methods used in [1] and [4] are also employed to 
compute the capacitance matrix of the same structure in [3]. According to [3], the methods 
used in [1] and [4] resulted in nonphysical values, for instance, negative self-capacitance 



C (nF/m) 


44 


£o 



w=l mm 

h— H 


Free space 


asms 


^ h=l mm 


Ground Plane 


Figure 3.2. A single microstrip line. 


Single Microstrip Case 



Figure 3.3. Comparison of the present method with the spectral domain approach as a 
function of the number of basis functions. 
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Figure 3.4. The plot of the surface charge density on the microstrip shown in Fig. 3.2. 
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Figure 3.5. Three conductors in a layered medium. All dimensions of conductors and 
spacings are identical. 
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Table 3.1. Comparison of the capacitance matrix for the three-conductor structure shown 
in Fig. 3.5 with [4]. 


Computation 


Delbare and Zutter [4] 


141.41 

-21.492 

-0.8952 


-21.492 

92.951 

-17.859 


-0.8952 

-17.859 

87.494 


(pF/m) 


142.09 

-21.733 

-0.8900 


-21.765 

93.529 

-18.097 


-0.8920 

-18.098 

87.962 


(pF/m) 


Potential Distribution 



Figure 3.6. The potential distribution of the system shown in Fig. 3.5 with the excitation at 
the center conductor. 
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Ground Plane 

Figure 3.7. Ten conductors in a layered medium. Dimensions of conductors are identical, 
and all units are in micrometers. 


Table 3.2. Capacitance matrix (pF/m) for the ten-conductor system shown in Fig. 3.7. 


307.13 

-41.2 

-11.34 

-6.28 

-5.351 

-218.8 

-4.966 

-1.385 

-0.814 

-0.729 

-41.20 

319.6 

-28.04 

-7.79 

-4.987 

-5.025 

-216.95 

-3.54 

-0.985 

-0.666 

-11.34 

-28.04 

310.5 

-24.29 

-8.587 

-1.366 

-3.503 

-218.4 

-3.18 

-1.148 

-6.279 

-7.79 

-24.29 

303.5 

-24.73 

-0.798 

-0.946 

-3.164 

-219.4 

-3.345 

-5.351 

-4.988 

-8.588 

-24.73 

290.5 

-0.708 

-0.627 

-1.12 

-3.325 

-221.3 

-218.8 

-5.01 

-1.363 

-0.796 

-0.709 

231.7 

-2.074 

-0.393 

-0.182 

-0.134 

-4.954 

-216.97 

-3.494 

-0.944 

-0.628 

-2.074 

232.0 

-1.19 

-0.255 

-0.135 

-1.382 

-3.532 

-218.4 

-3.157 

-1.12 

-0.393 

-1.19 

231.8- 

-0.8752 

-0.242 

-0.813 

-0.982 

-3.174 

-219.5 

-3.32 

-0.182 

-0.255 

-0.875 

231.6 

-0.763 

-0.729 

-0.665 

-1.145 

-3.34 

-221.4 

-0.134 

-0.135 

-0.242 

-0.763 

231.3 
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values. The capacitance values given in Table 3.2 are stable and agree well with the ones 
given in [3]. The method used in [4] required the CPU times of 8961 1.19 s on an IBM 
RS-6000 station with 300 basis functions, whereas the method in [3] took 458.67 s. On 
the same machine, the presented method took only 85.7 s of CPU time. 

For the final 2-D example, a multilayered stripline case, shown in Fig. 3.8, is also 
analyzed. Fifty basis functions are used to discretize each conductor. The computed 
capacitances were C n = C 22 = 217.07 pF/m and C 12 = C 21 = -107.76 pF/m. Data 
obtained from the Ansoft Maxwell software are C n = C 22 = 217.65 pF/m and C 12 = C 21 = 
-108.24 pF/m. 

3.3 Computation of the Inductance Matrix 

Inductances are obtained by solving magnetostatic problems, which, in general, are 
governed by a vector equation unlike electrostatic problems, which are governed by a scalar 
equation. As a consequence, the computations of inductances are often computationally 
more intensive than the computations of capacitances. However, for 2-D problems (the 
uniform transmission-line configurations) the magnetostatic problems can also be written 
by a scalar equation using only the z-component of the magnetic vector potential A . For 


Top Ground Plane 



Bottom Ground Plane 

Figure 3.8. Two conductors in a layered medium with two ground planes. Dimensions of 
conductors are identical. 
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these configurations, an isomorphism exists between the magnetostatic and electrostatic 
problems [1], [15], and the solutions of the magnetostatic problems are often obtained from 
solving the equivalent electrostatic problems. Therefore, the technique used to solve the 
capacitance in the previous section can still be applied to the inductance calculation. 

In the following subsection, this isomorphism between magnetostatic and 
electrostatic problems is first discussed, and then the expression relates the inductance 
matrix to the capacitance matrix of the equivalent electrostatic problem is derived. 

3.3.1 Isomorphism between electrostatic and magnetostatic problems 

To simplify the analysis, we will assume that there are two isotropic layers. Then, 
at the source-free region, electrostatic problems for a multiconductor system embedded in a 
multilayered medium are governed by the following Laplace equation: 

d^V d 2 V- ■? 

^- + ^- = V}Vi= 0 i = l,2 (3.8) 

dx 2 dy 2 ' 

where V t is the potential distribution in the ith layer, and the associated boundary condition 
at the dielectric interface is 


E{ = £ 2 A" = D 2 (3.9) 

Here, the superscripts t and n denote the transverse and normal components of the vectors. 
Rewriting the above equations in terms of V ( , 

(Z x n) • V, V, = (z x n) • V, V 2 n • (e rl V, V, ) = n ■ (e r2 V, V 2 ) (3.10) 

The associated boundary condition of (3.8) at the conductor interface is 

E \ — 0 Df = -q s (3.11) 

where q s is the surface charge density, and again, rewriting (3.1 1) in terms of V- , 

(£xn)-V,V* = 0 n(£„V,M) = -— (3-12) 

^ O 

Now for magnetostatic problems, we have the following equation at the source-free 


region: 
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V,xB,= 0 « = 1,2 (3.13) 

where fi, is the magnetic flux density in the ith layer. Using the facts, B, = V x A, and 
only the z-component of the magnetic vector potential A,- is nonzero under the quasi-static 
approximation, the above equation can be written as 

V, x (£ x VfAj,) = z(V, • V f A z/ ) = 0 => V^, = 0 (3. 14) 

The associated boundary condition of (3.13) at the dielectric interface is 

H{ = H2 = (3.15) 

and the above equations in terms of Aj, are 

n(— V f A zl ) = n-(— V,A z2 ) n(zxV f A zl ) = n(zxV,A z2 ) (3.16) 
Prl Pr2 

Applying vector identities to the above equations, we have 

n (— V f A zl ) = n-(— V^) (t xn)- = (? xn) V f ^ 2 (3.17) 
Mrl ^2 

The associated boundary condition of (3.13) at the conductor interface is 

nxHi = J z z B?= 0 (3.18) 

where J z is the longitudinal surface current density, and again, rewriting the above 
conditions in terms of A zb 

n(— V,A zt ) = -// 0 7 z (z xn)- V^Aj, =0 (3.19) 

Pri 

Now comparing (3.8), (3.10), and (3.12) with (3.14), (3.17), and (3.19), it is 
clear that if we replace l//r n and p 0 J z with e ri and q s /e 0 in (3.17) and (3.19), A zi is 
equal to V t . Therefore, A ZI can be computed by solving the equivalent electrostatic 
problem with replacing £ n by 1 / fi ri , and the current distribution J z (p) can be obtained 
using the following formula: 

J z (p) = c 2 q s (p) 


(3.20) 
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where c is the speed of light, and q s (p) is the surface charge density of the equivalent 
electrostatic problem. 

The inductance matrix L is defined to relate the magnetic flux differences between 
the signal conductors and the reference conductor to the currents on the signal conductors 
in the following manner: 

Li = \jr (3.21) 

where the ith elements of y/ and / are the magnetic flux difference between the ith signal 
conductor and the reference conductor and the current on the ith conductor. The inductance 
matrix can now be determined by solving the current distributions for N independent y/ t 
and can be put into the following matrix form: 

L = \|/I" 1 (3-22) 

where the ith columns of \\f and I are the ith magnetic flux excitation vector and current 
vector. Solving the equivalent electrostatic problem with V = \|f, and using the fact that 
I = c^Qeq , we have 

L = 4 C «q (3 ’ 23) 

c 

where Qeq and C eq are Q and C matrices in (3.4) for the equivalent electrostatic problem. 
For nonmagnetic media, the above equation becomes 

L-- Ic? (3.24) 

c 

where C 0 is the capacitance matrix for a free-space case in which all dielectric layers are 
replaced with free space. 

Since the inductance matrix is related to the inverse of the capacitance matrix, the 
following properties can be deduced from (3.5a), (3.5b), and (3.5c): 

Lij = Lji (3.25a) 

kj > 0 


(3.25b) 
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ki>^kj (3.25c) 

j*i 

Hence, L is also a positive definite symmetric matrix. This Maxwellian inductance matrix 
defined by (3.21) is related to the physical self- and mutual inductances by 

kj = kj (3.26a) 

li,i = L s (3.26b) 

3.3.2 Numerical examples 

Since the inductance matrix can be obtained using the electrostatic solution 
technique discussed in Section 3.2 and this technique has been shown to be very accurate 
and efficient, the resulting inductance matrix from this technique is also expected to be 
accurate and efficient. Hence, only a few examples are presented for the inductance 
calculations in this section, and more examples will be given as the remaining two 
parameters, the conductance and resistance matrices, are discussed. 

The three-conductor system shown in Fig. 3.9 is considered as a first example. 
The resulting capacitance and inductance matrices are compared with others in Table 3.3. 
The maximum number of exponential functions used in the approximation was seven, and 
16 basis functions are used to discretize each rectangular conductor, whereas 24 basis 
functions are used for the circular conductor. For a second example, the three-conductor 
system shown in Fig. 3.5 in Section 3.2 is again considered, and a comparison with the 
result obtained from [4] is given in Table 3.4. An excellent agreement was found for both 
cases. 

3.4 Computation of the Conductance Matrix 
3.4.1 Theory 

As mentioned in the beginning of this chapter, the inclusion of dielectric losses does 
not violate the quasi-static condition of a multiconductor system, and the conductance 
matrix G, which models losses due to the finite loss tangent of an imperfect dielectric, can 
be computed from the solution of the Laplace equation by means of a complex dielectric 
constant [10]. Since the dielectric constant is complex for this approach, the closed-form 
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Figure 3.9. Three conductors in a layered medium. The dimensions used in coordinates 
are in millimeters. 


Table 3.3. Comparison of the inductance matrix for the three-conductor structure shown 
in Fig. 3.9. 


C (pF/m) 


L (nH/m) 


Computation 


Wei et al. [1] 


Delbare and Zutter [4] 


116.9 

-13.33 

-62.79 

-13.33 

33.71 

74.96 

-62.79 

74.96 

362.3 

124.4 

-13.00 

-68.25* 

-13.00 

33.40 

-7.196 

-68.25 

-7.196 

352.3 _ 

125.9 

-13.12 

-69.55* 

-13.12 

34.10 

-7.182 

-69.55 

-7.182 

357.6 


489.0 

196.1 

114.7 

196.1 

612.8 

74.96 

114.7 

74.96 

219.1 

*496.5 

199.6 

118.3* 

199.6 

616.3 

77.28 

118.3 

77.28 

233. 1_ 

491.9 

198.9 

117.5 

198.9 

612.83 

76.781 

117.5 

76.781 

229.94 








Table 3.4. Comparison of the inductance matrix for the three-conductor structure shown 
in Fig. 3.5 in Section 3.2 with [4]. 


Computation 


Delbare and Zutter [4] 


277.90 

87.569 

36.629' 


'277.73 

87.758 

36.770' 

87.569 

328.87 

115.52 

(nH/m) 

87.758 

328.60 

115.77 

36.629 

115.52 

338. 16_ 


36.770 

115.77 

337.98_ 


(nH/m) 


Green’s function discussed in Chapter 2, which was obtained using the real exponential 
approximation, can no longer be applied. In principle, a new closed-form expression for 
the Green’s function which accounts for the complex dielectric constant can be derived; 
however, this approach will be computationally intensive as it involves complex 
operations, and a simple perturbation approach is employed in this section. 

The conductance matrix G for a multiconductor system embedded in a multilayered 
lossy dielectric medium is defined by 

G Vi = Ii i = l,2,- -,AT (3.27) 

where V-'s are N independent line voltage vectors. If' s are the corresponding shunt 
current vectors per unit length that arise due to dielectric losses, and N is the number of 
modes (N signal conductors or N lines). Equation (3.27) in the matrix form is 

G = I S V-' (3.28) 

where the ith columns of matrices V and I s are V* and If, respectively. The jith element 
of the matrix I, which is the shunt current on the jth line due to the voltage excitation of V* , 
can be obtained from the normal component of the electric field at the surface of the 
corresponding conductor of the jth line: 
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Here, Jf(p ) and £)(p) are the shunt surface current density and the electric field 
corresponding to the excitation vector, a(p) are the conductivities of the dielectric 
media, h is the surface normal vector, and Cj is the surface contour of the conductor. 
Using the fact that the normal component of the displacement field at the surface of a 
conductor is related to the charge density, (3.29) can be rewritten as 

ij,i = ig l ip' )<y(P' )/e(p' W (3.30) 

Jc j 

It can be shown that the potential distribution near the surface of the conductor does 
not change substantially when a slight loss is introduced in the dielectrics'; consequently, 
neither the electric field nor the charge density is much different for the lossy and lossless 
cases on the conductor. Thus, the charge density from a lossless system, which can be 
determined using the method given in Section 3.2, can still be used as the charge density 
q'(p) of a lossy system. Finally, once the charge densities have been obtained for N 
independent voltage excitations from the lossless system, the conductance matrix G can be 
derived from (3.28) and (3.30). 

Since the I matrix in (3.28) is directly related to the Q matrix in (3.8), the properties 
of the conductance matrix are the same as that for the capacitance matrix: 

Gij = G u (3.31a) 

Gjj > 0 and G, ; < 0 for j * i (3.31b) 



> 



(3.31c) 


and the relationship between the conductance matrix and the physical self- and mutual 
conductances is similar to that for the capacitance case, and is given by 

Gij=-G% (3.32a) 


'This can be seen from the fact that neither the lines of force nor the equipotential lines near the 
surface of the conductor are affected by the introduction of the dielectric losses. 
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G m =G ' + X G 5 (3.32b) 

j*i 

3.4.2 Numerical examples 

A cylindrical conductor above a perfectly conducting ground plane shown in Fig. 
3.10 is considered for a first example. Analytical formulas for the tr ansmis sion line 
parameters of this system can be obtained from the two-wire line case [1], [17] and given 
by 


:^-cosh _1 (2 Hid) 

2lt 

(3.33a) 

l7te 


cosh -1 (2H/ d) 

(3.33b) 

2R s H/d 

(3.33c) 

~m^j(2H/d) 2 -l 

2ko 


cosh -1 (2///d) 

(3.33d) 


where H is the distance from the center of the cylinder to the ground plane, d is the 
diameter of the cylinder, R s is the surface resistivity of the conductor, and cr is the 
conductivity of the dielectric medium. The computed results are compared in Table 3.5 
with the analytical solutions and numerical results from [1]. Forty basis functions were 
used to solve the electrostatic problem. The comparison shows excellent agreement. 

For the next example, the four-conductor system shown in Fig. 3.1 1 is considered, 
and the results are compared in Table 3.6. The loss tangents of the bottom, middle, and 
top dielectric layers are 6.4 x 10* 4 , 8.6 x 1(H, and 2.4 x 1(H, respectively. 
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Figure 3.10. A cylindrical conductor above a PEC ground plane. 


Table 3.5. Comparison of the computed data with the analytical solutions for the 
cylindrical case shown in Fig. 3.10. 

Analytical Solutions Computed Data Wei et al. [1] 

C (pF/m) 89.80 89.75 89.44 


G (pS/m) 


67.71 


67.66 


67.44 
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Table 3.6. Comparison of the computed results with data obtained from [1] for the four- 
conductor system shown in Fig. 3.11. 
















3.5 Computation of the Resistance Matrix 
3.5.1 Theory 
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At the frequencies near dc, the current is uniformly distributed throughout the cross 
sections of conductors, and the resistance matrix, which is diagonal, is ob tain ed by takin g 
the inverse of the product of the conductivity and the cross-sectional area of each 
conductor. In this section, we will focus on the computation of the resistance matrix in the 
high-frequency region where the skin effect is prominent, and the transverse current 
distribution is approximately the same as that for a perfect conductor case. 

At this high-frequency region, the transverse current distribution on lossy 
conductors is no longer uniform due to the edge and the proximity effects [16], [17], and 
the longitudinal voltage drop due to the surface resistivities of the conductors becomes a 
function of position in the transverse plane. This nonuniform voltage drop cannot be 
characterized by the resistance matrix; hence, the effect of the finite conductivities of 
conductors is rather complicated in the high-frequency region in contrast to the effect of the 
finite loss tangents of the dielectric media, which was readily modeled with the conductance 
matrix G based on the shunt current. However, if one is only interested in modeling the 
power loss on conductors, the resistance matrix R can still be used to characterize the 
conductor power loss per unit length, and is computed either from the conductor power 
loss or from attenuation constants, which, in turn, are obtained from the power loss on 
conductors [10], [11]. Modeling the power loss with the resistance matrix results in the 
effective longitudinal voltage drop, which gives the correct power loss per unit length and 
is uniform over the transverse plane 

Unfortunately, the power loss on the conductor is not only dependent on the 
magnitude of the current but is also a function of the current distribution on the surface of 
the conductor; therefore, the resistance matrix is also affected by the cuirent distributions 
employed to compute it. Furthermore, although the current distribution for a single 
transmission line is not affected by excitations except for scaling of the magnitude, it 
becomes strongly dependent on the types of excitations for coupled transmission lines [11], 
Hence, it follows that the resistance matrix, in general, is dependent on excitations for 
multiconductor systems. Because of this undesirable dependency, the resistance matrix 
may not be computed unless the current distribution has been determined. In this section, 
the nondiagonal resistance matrix is first considered to illustrate this dependency; then, the 
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diagonal resistance matrix, which is relatively insensitive to the choice of excitations, is 
proposed. 


To simplify the analysis, let us assume that there is only one ground plane at the 
bottom of the dielectric layers and that the ground plane is perfectly conducting. Then, let 
us define the resistance matrix Rj such that it gives the correct power losses on each 
conductor as follows: 

Rj if = yf i = 1,2,- • -,N (3.34) 


where if' s are N independent conduction current vectors, the y'th element of if is given 
by 



(3.35) 


and Vffi' s are the line voltage vectors, which represent the effective voltage drops, and are 
computed from the power loss due to the conduction current vectors if' s . The jth 
element of Vf^ is 






(3.36) 


where Pij is a power loss on the ;th conductor, and R s j is the surface resistivity [17] of 
the jth conductor given by 


^s,j — . 


w 


°J s j 


(3.37) 


Here, S is the skin depth. The above equation is valid when the thickness of a conductor is 
at least two skin depths thick, and it is assumed in this section that (3.37) is valid for the 
frequency of interest. 

Now the resistance matrix can be readily obtained once the conduction surface 
current densities Jf(p), which result in N independent conduction current vectors, are 
determined. In the high-frequency region where the skin effect is prominent, the surface 
current density Jf(p) on good conductors is close to that for perfect conductors. 
Moreover, as discussed in Section 3.3, the surface current density on perfect conductors 
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can be obtained by solving the equivalent electrostatic problem, which is constructed by 
replacing all dielectric media with free space. Hence, Jf(p) can be computed using 
(3.20). This approximation of Jf (p) with the current density for a perfect conductor case 
is also valid for cases where losses are not relatively large but they are mostly due to the 
small cross sections of conductors instead of the small values of finite conductivities [18]. 

Equation (3.34) is applied to compute the resistance matrix for the coupled 
microstrip lines shown in Fig. 3.12. The conductivity of both microstrip lines is 
57.6 MS/m, and the ground plane is assumed to be perfectly conducting. The method 
described in Section 3.2 is used to solve the equivalent electrostatic problem, and forty 
pulse basis functions per each line are used to represent the charge distribution. The 
computed values of Rj at a frequency of 1 MHz with various excitations are shown in 
Table 3.7. The first column shows the excitation voltage used in the equivalent electrostatic 
problem, and the second column shows the resulting excitation currents. The third to sixth 
columns correspond to the computed values of Rj based on (3.34). All elements of Rj 
varied as the excitation vectors are changed except for the case in which the two excitation 
vectors are simply scaled (see first and second rows); the variation is especially prominent 
for the off-diagonal elements of Rj, which changed from negative to positive when 
different excitations were used. 

As observed from Table 3.7, the diagonal elements of the resistance matrix are 
somewhat less sensitive to the choice of current excitations than the off-diagonal elements. 
Hence, it may be natural to define the resistance matrix as a diagonal matrix, as it is for a dc 
case, such that the resistance matrix becomes less dependent on the choice of excitations. 
To obtain the diagonal resistance matrix, let us define the resistance matrix by matching the 
total power loss on the multiconductor system rather than by matching the power losses on 
each conductor: 


l[VL 2 li= p i i = 1,2,- • -,N (3.38) 

T 

where R 2 is assumed to be a diagonal matrix, j is the transpose of if, and P, is the 
total power loss on the system due to the current vector if, obtained by 
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Figure 3.12. Coupled microstrip lines. 


Table 3.7. The resistance matrix R of the coupled microstrip lines shown in Fig. 3.12 
with various excitations. 


Excitation 

R (nondiagonal) (|iQ/m) 

R (diagonal) 
((lQ/m) 

Voltage 

Current 

Rn 

Rl2 

R21 

R 22 

R 11 

R22 

(1,0) 
(0, 1) 

(2.00e9, -7.2 le8) 
(-7.2 le8, 2.00e9) 

23.65 

-5.439 

-5.439 

23.65 

27.13 

27.13 

^^9 

(9.99e9, -3.60e9) 
(-1.44e9, 4.00e9) 

23.65 

-5.439 

-5.439 

23.65 

27.13 

27.13 

(i, i) 
(i,-D 

(1.28e9, 1.28e9) 
(2.72e9, -2.72e9) 

26.38 

-1.176 

-1.176 

26.38 

• 

• 

(i, i) 
(1,2) 

(1.28e9, 1.28e9) 
(5.56e8, 3.27e9) 

17.00 

8.200 

0.416 

24.78 

24.37 

26.03 

| 

(2.00e9, -7.21e8) 
(5.56e8, 3.27e9) 

27.91 

6.347 

-4.717 

25.66 

27.28 

25.95 















where N is the number of signal conductors in the system. Now, noting that R 2 is a 
diagonal matrix, (3.38) can be put into the following matrix form: 

Z R = P (3.40) 

Here, the ith element of the vector P is Pi, and R is a vector which consists of the 
diagonal elements of R 2 . The i/th element of Z is a square of if j, where if j is the y'th 
component of If . The diagonal resistance matrix can be obtained by solving the above 
system of linear equations. Note that N independent conduction current vectors of if no 
longer guarantee the matrix Z to be nonsingular. One simple way to obtain N current 
vectors, which result in a nonsingular Z matrix, is to excite the voltage only at one 
conductor at a time in the equivalent electrostatic problem as shown in the first row in Table 
3.7. This choice of excitation vectors is shown to result in a more stable resistance matrix 
as shown in the next numerical example. 

To examine the variation of the diagonal resistance matrix for the different choices 
of excitation vectors, the previous microstrip lines are considered again, and the computed 
results are shown in Table 3.7. As discussed in the previous paragraph, when the voltage 
excitations of (1, 1) and (1, -1) are used, Z becomes a singular matrix; the diagonal 
resistance matrix could not be determined for this case. As the resistance matrix is defined 
to be a diagonal matrix, the nondiagonal elements, which are strongly affected by the 
choice of excitation vectors, no longer cause any problem in the diagonal resistance matrix. 
Furthermore, as shown in Table 3.7, the diagonal elements of the diagonal resistance 
matrix are also less dependent on the choice of excitation vectors than the diagonal elements 
of the nondiagonal resistance matrix. The best choice of excitation vectors appears to be 
the first two rows in Table 3.7 since the two elements of the resistance matrix are identical 
as expected from the symmetry of the system. 

In practice, the excitations of the resistance matrix based on (3.34) should be 
chosen so that the resulting resistance matrix is a symmetric matrix to preserve the 
reciprocity principle. Being a diagonal matrix, reciprocity is naturally satisfied for the 
resistance matrix based on (3.38). In general, the Maxwellian resistance matrix calculated 
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from (3.34), (3.38), or the perturbational analysis on attenuation constants is related to the 
physical self- and mutual resistances by 

RiJ = Ri!j (3.4 i a) 

R U =R S (3.41b) 

3.5.2. Losses due to imperfectly conducting ground planes 

To account for losses due to imperfectly conducting ground planes, the surface 
current distribution on the ground planes must be determined for the given current 
excitation on the signal lines. In general, an integral equation has to be formulated to 
compute the surface current density on the ground planes, and a detailed discussion of the 
determination of the current densities on ground planes is presented in this section. 

Let us first assume that a multilayered dielectric medium is backed by two ground 
planes below and above the layers, as shown in Fig. 3.13. Let J(x,y) , JS- ] (x ) , and J&’ 2 (x) 
be the current densities on the surfaces of the signal lines, the bottom ground plane, and the 
top ground plane, respectively. Then, assume that the surface current density on the signal 
lines J(x,y) is given by 

N, 

J(x,y)=^J k P k (x,y) (3.42) 

fc=l 


where P k (x,y) is an appropriate pulse basis function used to expand the current density on 
the surface of signal conductors, and N t is the total number of basis functions used to 
discretize the contours of the conductors. Then, the tangential magnetic field H^x) at the 
bottom ground plane due to J(x,y) and JS> 2 (x) at the bottom ground plane can be written as 
(see Fig. 3.13) 


1 N, 

#,(*>=— y 




2 nf^(x-x k ) 2 + y k 2x J ~(x-x’) +h 


-f 

7ir J- 


J g2 (x' )h 


-(be 


(3.43) 


Since the tangential magnetic field due to all the image currents of J(x,y) and JS> 2 (x) 
must also be equal to Hrfx) at the bottom ground plane, the total tangential magnetic field at 
the bottom ground plane must be twice that of H^x ) . Since the negative value of this total 
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Figure 3.13. The fcth segment of conductor contours between the two ground planes. 


tangential magnetic field is equal to the actual surface current density JS ] (x), the following 
integral equation can be formulated: 


j*Hx) + -f 00 j8 ' 2{x 2 )h 2 dx' = -I y — htyk 

n J ^( x -x') 2 +h 2 K ^x-x k ) 2 + yl 


(3.44) 


A similar argument holds for the tangential magnetic field at the top ground plane, and the 
resulting integral equation is 



Jg '\x')h _ ly WltZlkl 
(x-x') 2 +h 2 n£-‘( x -x k ) 2 +(h-y k ) 2 


(3.45) 


Now, truncating the ground planes to finite lengths and expanding J8- ] (x) and J8< 2 (x) 
using the pulse basis functions, 2 (3.44) and (3.45) can be put into the following form: 


I 

M 2 i 


m 12 

I 



(3.46) 


2 This truncation is possible since the actual current is confined near the conductor regions. To 
achieve good computational efficiency, nonuniform pulse basis functions should be used in the expansion 
of the current. 
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where Mjj is an N gi by Ngj matrix, JB»' and Vj are N gi dimensional vectors, N gi is the 
number of basis functions used to expand JS>‘(x), and I denotes the identity matrix. The 
expressions for the elements of the matrices and vectors are given by 



1 ifh 

Jt (xf 1 -x 8 j' 2 ) 2 +h 2 




ifh 




_ 


xf) 2 +h 2 


(3.47) 




hkyk 

- 4) 2 


+yl 



Jkhib-y k ) 

- x s k ) 2 +(h-y k ) 2 


(3.48) 


where l?’ m is the length associated with Jf m , which is the jth element of J8’ m , and x m 
and y m represent the center point of the mth basis function. The superscripts g, 7, g,2, and 5 
are attached to Jt to denote x’s on the bottom and top ground planes and the signal 
conductors, respectively. Now performing a simple substitution, a computationally more 
efficient form of (3.46) can be obtained: 

[I - M 21 M 12 ]J 8 ’ 2 = V 2 - M 21 Vi (3.49) 

J8’ 1 =Vj — M 12 J g ’ 2 (3.50) 


Now, the current distribution on the ground planes can be obtained by first solving the 
system of linear equations given by (3.49) for J8* 2 and using (3.50) for J& 1 . 

When there is no top ground plane, (3.46) is simplified to 

j 8,1 = Vj (3.51) 


Note that no solution of a system of linear equations is necessary for this case. 

Once the current density on the ground planes has been determined, (3.39) can be 
modified to include ground losses as follows: 

N N g 2 

Rs,j[ J i(P')f d P' + ^§ R s,{ J i'4p')\dp' (3.52) 

j=\ Jc J j - 1 Cj 

where Ng is the number of ground planes, which is either one or two. 
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3.5.3 Remarks on the various methods used in the computation of the 
resistance matrix 

In this section, several remarks on the various methods used to compute the 
resistance matrix are given. Since power losses are nonlinearly dependent on the current 
distribution for a multiline case, the resistance matrix becomes dependent on the 
excitations. Thus, if power losses are modeled with the resistance matrix, the resistance 
matrix cannot be independent of the excitations used in the computation. For a single-line 
case, the resistance matrices based on (3.34), (3.38) and attenuation constants are all 
identical, and they model the exact power loss per unit length regardless of excitations. For 
a multiline case, the nondiagonal resistance matrix defined by (3.34) is only valid for the 
case in which the actual current distribution belongs to one of the current excitations used in 
the computation. Similarly, the conventional approach based on the perturbation on 
attenuation constants [10] uses the modal power losses; therefore, the resulting resistance 
matrix is only valid for the modal current. However, the currents on the conductors are 
arbitrary so any current distribution is as important as the modal current distribution. 
Hence, the nondiagonal matrices based on (3.34) and the conventional approach are 
inappropriate for the computation of the resistance matrix for a multiconductor case. The 
diagonal resistance matrix based on (3.38) is shown to be least affected by the choice of 
current excitations and, hence, is most suitable for the computation of the resistance matrix. 

Now let us discuss the accuracy of the various definitions of the resistance matrix. 
First, the nondiagonal resistance matrix based on (3.34), which, in fact, is included in the 
chapter only to demonstrate the nonlinear effect, is symmetric for the modal current 
excitations (the third row in Table 3.7) and more accurate than the one from [10] since it 
matches the power losses on each conductor instead of matching the total power loss. 
Thus, (3.34) with the modal excitations can be used to define the unique resistance matrix, 
and this definition of the resistance matrix is already more accurate than the conventional 
approach [10]. The diagonal resistance matrix based on (3.38) also models the total power 
losses for chosen current excitations; therefore, it is as accurate as the one from the 
conventional approach. Conclusively, if the excitation is fixed as the conventional 
approach, the resistance matrices based on (3.34) and (3.38) are as accurate as the one 
from the conventional approach and unique. Although the resistance matrix based on 
(3.38) is less accurate in terms of modeling of power losses than the one based on (3.34) 
for a fixed excitation, it is more accurate for an arbitrary current distribution because it 
varies less for different excitations. Similarly, the resistance matrix obtained from the 
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conventional approach is also nondiagonal and, hence, is less accurate than the diagonal 
resistance matrix for an arbitrary current distribution. 

3.5.4. Numerical examples 

In this section, the diagonal resistance matrix for various geometries are computed. 
Most examples are extracted from [10] for comparison purposes. The thicknesses of the 
ground planes are assumed to be infinite for all examples. The method described in Section 
3.2 is used to solve all electrostatic problems. As a first example, the resistance value of 
the cylindrical conductor above a PEC ground plane shown in Fig. 3.10 on p. 55 is 
computed. The computed value was 0.8459 ft/m at a frequency of 100 MHz, whereas the 
closed-form formula (3.33c) was 0.8426 ft/m. Forty pulse basis functions were used to 
solve the electrostatic problem in our approach. As the next example, a single microstrip 
line is considered. The same geometry used in Section 3.5.1 (see Fig. 3.12) is analyzed 
with only one line. At a frequency of 1 MHz, the computed resistance value was 
2.493 x 10' 5 ft/m without any ground loss and 2.688 x 10- 5 ft/m with ground loss. The 
ground plane is assumed to be of the same material as the microstrip line, both materials 
having a = 5.76 x 10 7 S/m. The resistance value from [10] with ground loss was 
2.400 x IO- 5 ft/m after applying the proper scaling factor. The current distribution on the 
ground plane is shown in Fig. 3.14. The resistance values were also computed for a wide 
range of frequencies, and the results are plotted in Fig. 3.15. 

For the final example, the four-conductor system shown in Fig. 3.11 on p. 56 is 
considered, and the results are compared in Table 3.8. The current densities on the top and 
bottom ground planes with the excitation on the first conductor are shown in Figs. 3. 16(a) 
and (b). The resistance matrix shown in Table 3.8 is actually incorrect since (3.37) for the 
surface resistivity R s is no longer valid because the thicknesses of the first and third 
conductors are smaller than the skin depth 8. To obtain an approximate solution, we have 
defined the effective depth d e ff to be either 8 or the half thickness of the conductor, 
whichever is smaller, and the surface resistivity R s is redefined as 

1 

Odeffip) 


R S (P) = 


(3.53) 



































Jz (A/mm) Jz (A/mm) 
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The Current Density on the Top Ground Plane 



x (mm) 
(a) 


The Current Density on the Bottom Ground Plane 



x (mm) 
(b) 


Figure 3.16. The current densities on (a) the top and (b) the bottom ground planes for four- 
conductor system with the excitation on the first conductor. 
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The half thickness of the conductor is used since the surface current exists on both sides of 
the conductor. Note that the effective thickness d e ff is a function of p since the thickness 
of a conductor, in general, changes as p varies. The resistance matrix R for the previous 
case is computed again using the above formula, and the following result shows that the 
resistance value is substantially increased: 


R (|ifi/m) = 


4688 

0 

0 

0 


0 

836.9 

0 

0 


0 

0 

2539 

0 


0 

0 

0 

598.4 


3.6. Summary 

Accurate and efficient ways to compute the transmission line parameters of a 
multiconductor system embedded in a multilayered dielectric medium were presented in this 
chapter. The capacitance matrix was computed based on the closed-form Green’s function 
discussed in the previous chapter; thus, no numerical integrations or nested infinite 
summations are involved in the computation. The inductance matrix is calculated by 
solving the equivalent capacitance problem. The solutions (charge distributions) obtained 
from the computations of the capacitance matrix and the inductance matrix were directly 
used to compute the resistance and conductance matrices, rather than by using the usual 
perturbation analysis on the modal power. Hence, no additional electrostatic problem was 
solved for the computations of the resistance and conductance matrices. The diagonal 
resistance matrix was proposed in this chapter; this diagonal matrix has been shown to be 
relatively insensitive to the choice of the current excitations as opposed to the sensitivity of 
the traditional nondiagonal resistance matrix. Hence, it is more suitable for the computation 
of the resistance matrix. In addition to losses on signal conductors, those due to 
imperfectly conducting ground planes were also incorporated into the resistance matrix. 
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CHAPTER 4 

COMPUTATION OF THE EQUIVALENT CAPACITANCES OF VARIOUS 
STRIP DISCONTINUITIES: AN OPEN END, A BEND, 

AND VARIOUS JUNCTIONS 

4 . 1 Introduction 

Quasi-static analysis is often performed to characterize strip discontinuities, such as 
an open end, a bend, a T junction, and a cross junction, when the dimensions of the 
discontinuities are much smaller than the wavelength. A brief introduction of the quasi- 
static approximation was presented in Section 1 .2. Under this quasi-static analysis, the 
dominant effect of strip discontinuities is fringing fields due to the physical irregularities of 
discontinuity geometries, and the modeling of these fringing fields in terms of an equivalent 
(excess) capacitance is discussed in this chapter. 

Numerous papers have been published to compute the excess capacitances of 
various microstrip discontinuities, and a summary of popular methods can be found in [ 1 ], 
[2]. The usual approach involves finding the charge distributions in the presence of a 
discontinuity and in the absence of a discontinuity. Total charges are then obtained from 
these charge distributions, and the excess capacitance is computed by subtracting these total 
charges. This approach is referred to as the total charge formulation in this thesis, and it is 
widely known to be inaccurate due to the numerical error associated with the subtraction of 
the total charges due to reasons explained in [l]-[3]. Since the discontinuity effects are 
rather localized, the total excess charge distribution of the discontinuity system is often 
much smaller than the total charge. Hence, the total charges of the system with and without 
the discontinuity are nearly equal, and the errors associated with these two total charges are 
relatively very large compared to the total excess charge, which results in the equivalent 
capacitance. Moreover, the total charge formulation is also computationally expensive for 
the following two reasons. First, two problems must be solved for one with a 
discontinuity and one without a discontinuity. Second, discontinuities are associated with 
semi-infinite lines, and after truncating these lines to finite lengths, the unknowns 
associated with the total charge distribution must be placed over the whole surface of the 
truncated lines although the excess charge is localized around the discontinuities. 

The formulation of an integral equation in terms of the excess charge distribution, 
first proposed by Silvester and Benedek [3], has been applied to analyze various microstrip 
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discontinuities [3]-[6]. This approach overcomes the previously discussed numerical 
problems associated with the total charge formulation since the excess charge distribution is 
directly modeled as an unknown in an integral equation; in other words, the accuracy 
resulting from solving an integral equation is preserved for the final answer. This approach 
is referred as the excess charge formulation throughout this thesis. The Green’s function 
for a layered medium is employed in this approach. For N dielectric layers, the expression 
for this Green’s function would consist of an N-l nested infinite series [7]; hence, in 
practice, this form of the Green’s function may not be applied to a multilayered medium 1 . 
Recently, Sarkar et al. [8] solved the discontinuity problems for a multilayered medium 
using the free-space Green’s function, but additional unknown charges (over unknown 
charges on the surface of a conductor) had to be placed on the dielectric interfaces and the 
top ground plane to model the polarization charge and the free charge. Although the 
inclusion of these additional unknowns may be tolerable for 2-D problems, it is 
computationally too burdensome for 3-D problems. 

In this chapter, the closed-form Green’s function discussed in Chapter 2 is 
employed to formulate an integral equation in terms of the excess charge distribution; thus, 
the presented method requires neither additional unknowns to model dielectric interfaces 
and the top ground plane nor evaluations of any infinite series except for cases where the 
top ground plane is present. When the top ground plane is present, using the closed-form 
Green’s function is still numerically advantageous since the nested infinite series in the 
expression of the usual Green’s function become a simple infinite series without nesting. 
Although it is possible to avoid infinite series even for cases where the top ground plane is 
present by modeling it as an additional conductor, 2 it substantially increases the number of 
unknowns; hence, it is not considered in this chapter. 

In Section 4.2, the general description of discontinuity structures is discussed, and 
the general representation of strip discontinuities, from which most of common 
discontinuity types can be derived, is presented. In Section 4.3, the closed-form Green’s 
function is employed to formulate an integral equation to determine the excess charge 
distribution of the generalized discontinuity structure given in Section 4.2; thus, only one 


'See also Section 2.2 for more detailed discussion of this Green s function. 

2 This approach is used in Chapter 6 for the computation of the equivalent capacitance of a strip 
crossover. 
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integral equation is formulated for various types of discontinuities, unlike other approaches 
given in [3]-[6], in which a different integral equation has been formulated for each 
discontinuity type. Then, the method of moments is employed to solve the integral 
equation. The equivalent capacitances of an open end, a step junction, a bend, and a T 
junction are considered as numerical examples, and the computed capacitances are 
compared with other published results in Section 4.4. 

4.2 General Statement of the Problem 

The planar view of the general geometry of a discontinuity considered in this 
chapter is shown in Fig. 4.1. A discontinuity, in general, consists of an arbitrary number 
N t of traces, which are all connected either with or without a junction region. This general 
geometry represents most of the common strip discontinuities, e.g., an open end, a 
nonorthogonal bend, and various junctions. Although the present approach can handle 
conductors with finite thicknesses, the conductor thicknesses are assumed to be infinitely 
thin in this chapter. 

The discontinuity structure is embedded in a layered dielectric medium, which is 
shown in Fig. 4.2. An arbitrary number of dielectric layers are located on top of a ground 
plane, and the layered dielectric medium is terminated by an optional ground plane on the 
other side. Unlike the dielectric medium considered in Chapter 2 (see Fig. 2.1), it is 
assumed that the bottom ground plane always exsists. 

Some common strip discontinuities which can be represented with the general 
discontinuity geometry shown in Fig. 4. 1 are depicted in Fig. 4.3 with their equivalent 
circuit representations. These representations include transmission lines; hence, they 
assume that there is at least one ground plane. In general, the equivalent inductances 
should be added in these equivalent circuit representations, and they can be computed using 
the method described in [7]. 

4.3 Formulation of an Integral Equation 

The integral equation relating the electrostatic potential 0(r) and the charge density 
q{r) on the surface of a conductor for 3-D problems is given by 
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Figure 4. 1 . General geometry of a strip discontinuity. 


Optional Top Ground Plane 


y =d Nd 



Bottom Ground Plane 


x 


Figure 4.2. Cross-sectional view of a multilayered dielectric medium. 
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(d) A T junction 
Figure 4.3. Continued. 


<t>(r) = J G W (r\r' )q{r' )dr' = (g 3D ,q) (4.1) 

n 

where the i2’s are the surfaces of conductors: traces and a junction region. G (r|r' ) is 
the 3-D closed-form Green’s function for a layered medium, which accounts for the 
polarization charges and free charge on ground planes. To simplify the notation the 
integration is symbolically written as (•,•). Now let us rewrite the charge density q(r) in 
the following manner: 

f qj (r), if r is on the junction region 

q(r) = j . (4-2) 

[4r( r )’ if r is on the ith trace 

and, for each trace, let us decompose the charge densities q l j(r) into the uniform charge 
density q^'\r) and the excess charge density qj^ ess ' 1 {r) \ 

,>.(r)-^''(r)+ 9 f““ JJ (r) (4.3) 

Here, the uniform charge density q^’\r) is obtained by solving a 2-D problem, in which 
it is assumed that only the ith trace is present in the medium and that the ith trace is 
infinitely long in both directions. A detailed discussion for solving 2-D problems was 
given in Section 3.2.1. The uniform charge density qj nif '\r) exists only on the ith trace, 
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which is a semi-infinite line; hence, G semi (r\r 0 ,£) should be used to compute the potential 
due to qf*\r). 

Using (4.2), the integral equation (4.1) can be written as follows: 

<Kr) = (g 3Z> ,?) = (G^ gj ) + ^( G ^ T ) (4.4) 

i=\ 

Now using (4.3), the above equation is rewritten as 

^ N > 

0(r)-£{G“™^ ^ ^ r •V{c 3D , ^ ,) + y(c 3D .4«"«) (4.5) 

'=i ;=i 

A complete list of expressions of the closed-form Green s functions for a point charge, a 
line charge, and a semi-infinite line charge with or without a top ground plane is given in 
Section 2.3. It should be noted that all quantities in the left-hand side of (4.5) are known 
assuming 2-D problems have been solved a priori. 

Now applying the method of collocation with the excitation voltage <p 0 , the above 
integral equation yields the following system of linear equations: 


V 


‘ Vj 

' q“ nif ’ 1 ' 


' M ltl 

M l.N t +l 

_ excess, 1 
q T 

A. 


> t+ i_ 

q untf, N t 


_ M N t +l,l 

' M N t +l,N,+l_ 

excess, N, 

L^t J 


where 


[ M u]„,* = JJg 3 ^ V W (4.7a) 

[Vi ] p , q -jG semi (r p y4)dr 


(4.7b) 
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excess, i _ excess, i excess, i 

Hj “ HT % 1 ’*# 7\2 * 


excess, i 


(4.7c) 






unif,i _ Lum/,/ „unif,i 


q T - <? r i ,q T2 


unif ,i 
' q T.N 2D ^ 


(4.7d) 


(4.7e) 


Here, and C q denote the source patch and line segment for 3-D and 2-D problems, and 
r p is the observation (testing) point, which is the center of a patch or a line segment. Nsdj 
and Ny are the total numbers of patches to represent q p xcess ' l (r) and qj(r), respectively, 
whereas is the total number of line segments to represent qj' l {r). It is assumed 
that <(>o is also used as the excitation voltage in 2-D problems. The closed-form integration 
formulas for (4.7a) and (4.7b) were discussed in Section 2.4. Now given the excitation 
voltage 0 O , the excess charge distribution can be determined by solving the above linear 
system of equations. 


Now once (4.6) is solved, the excess (equivalent) capacitance Cf can be obtained 
by 


C e 


Nr 




i=i 



(4.8) 


where 


Nj 

Qj = '£‘U,k Area J,k 
*=i 


n j 


Q^^q^Area^ 


k = 1 


(4.9a) 


(4.9b) 


Here, Area y,* and Area l T k are the areas of the fcth patches used to represent qj(r) and 

qtxcess, i( r ), respectively. 

Throughout the formulation, we have assumed that the junction region exists 
between traces. The formulation for cases without the junction region, such as an open end 
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and step junctions, can be easily obtained by following the above procedure; in fact, (4.6) 
can be used by removing terms corresponding to qj{r). 

4.4 Numerical Examples 

A computer program is written based on the method discussed in the previous 
sections. The program assumes that the shape of the junction region (see Fig. 4.1), if it 
exists, is polygonal, and it can handle an arbitrary number of dielectric layers as well as 
traces. Each trace is represented by the center point of the ending side and the angle from 
the x-axis to the center line of a trace (see Fig. 4.4). Nonuniform meshing is utilized to 
reduce the number of unknowns. 

Excess capacitances for four common strip discontinuities, an open end, a step 
junction, a bend, and a T junction, shown in Fig. 4.3, are computed using the program. 
The following parameters are used: 1) an open end: w = 0.5tnm; 2) a step junction: 
wj =0.1 mm and w 2 =0.2 mm; 3) aright-angle bend: w l = w 2 =0.15 mm; and 4) a T 
junction: w ] =w 2 = w 3 =0.15 mm . Three different types of media are considered for each 
discontinuity with the following parameters (see Fig. 4.5): 1) an open end: e, =4.2, 
^2 — 2.5, yj =1.0 mm, y 2 =1.5 mm, and y 3 = 2.0 mm; 2) a step junction: £j = 6.0, 



Figure 4.4. Representation used to describe traces. 
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£ 2 =4.2, ^=0.1 mm, y 2 =0.2 mm, and >>3 =0.3 mm; 3) a bend: £j = 2.5, £ 2 =4.2, 
yj = 0.15 mm, y 2 = 0.3 mm, and >3= 0.5 mm; 4) a T junction: £j=2.5, £ 2 =4.2, 
yi = 0.15 mm, y 2 = 0.3 mm, and = 0.5 mm. All discontinuities are assumed to be 
embedded at y = yj . To place 3-D unknowns for the excess charge distribution, the length 
of each trace is truncated at l = 8w . The total numbers of unknowns per each trace were 50 
for a 2-D problem and 160 for a 3-D problem, whereas 100 unknowns were used for the 
junction region. The maximum number of exponential functions used to approximate each 
coefficient function Kf{y,m,n ) was 5. 

The computed results are shown in Table 4.1 with the comparison data for the 
microstrip case (Fig. 4.5(a)). A good agreement was found overall as shown in the table. 
It is interesting to note that for some cases the value of an excess capacitance turns out to be 
negative. Although a physical capacitance must be positive, an excess (equivalent) 
capacitance is hypothetical and can be negative. The excess charge distribution is plotted 
for the microstrip case in Figs. 4.6, 4.7, 4.8, and 4.9. The excess capacitance of an open 
end is also computed as a function of a trace width and compared in Fig. 4.10 with the 
closed-form formula from [2]. 

Since the number of unknowns can be quite large for 3-D problems, (4.6) is solved 
using the iterative methods, such as the Generalized Conjugate Residue (GCR) routine or 
the Generalized Minimum Residual (GMRES) routine instead of the LU factorization. The 
comparison of the CPU time used in the above computation for the microstrip case is given 
in Table 4.2. The CPU time is measured using an Alpha workstation, and the iteration is 
performed until the maximum relative error is smaller than 0.7%. The gain in the computer 
time was quite noticeable for the iterative approaches as the number of unknowns 
increases. Since (4.6) is diagonally dominant, both iterative approaches converged quickly 
as shown in Table 4.2. In general, the GCR routine performed better than the other 
appoaches; furthermore, it requires less computer memory than for the GMRES routine. 

4.5 Summary 

An efficient method to compute excess capacitances of strip discontinuities 
embedded in a multilayered dielectric medium was discussed in this chapter using the 
excess charge formulation in conjunction with the closed-form Green s functions; thus, the 
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Figure 4.5. Three media considered for numerical examples. 
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Table 4.1 . The numerical results (units are in femtofarads). 



Medium 1 

Medium 2 

Medium 3 

Computation 

Others 

Computation 

Computation 

Open End 

17.33 

17.0 [1] 

23.52 

19.62 

Step Junction 

1.120 

1.05 [1], 0.74 [8] 

1.352 

0.609 

Bend 

6.210 

6.75 [1], 5.8 [8] 

7.006 

9.184 

T Junction 

1.385 

1.9 [8] 

-4.917 

-0.818 



0.5 


l(mm) 


Figure 4.6. The excess charge distribution for the open end case. 
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Figure 4.7. The excess charge distribution on (a) the wider trace and (b) the narrow trace 
for the step junction case. 
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Figure 4.9. The excess charge distribution on (a) the side traces, (b) the center trace, and 
(c) the junction region for the T junction case. 
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Table 4.2. The CPU time used in the computations. 



GCR 

GMRES 

LU 

Time (sec) 

Iterations 

Time (sec) 

Iterations 

Time (sec) 

Open End 

12.42 

17 

12.54 

19 

12.71 

Step Junction 

73.40 

17 

73.49 

20 

76.37 

Bend 

88.39 

10 

8.67 

15 

95.81 

T Junction 

116.73 

9 

167.54 

16 

187.16 


presented method is accurate and numerically efficient. Unlike other approaches, only one 
integral equation was employed in this chapter to handle various strip discontinuities 
instead of formulating a different integral equation for each discontinuity type. The 
numerical results for the microstrip case agreed well with other published results. 
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CHAPTER 5 

COMPUTATION OF THE EQUIVALENT CAPACITANCE OF A VIA IN A 

MULTILAYERED BOARD 

5 . 1 Introduction 

Although a via is one of the most common discontinuities encountered in high- 
speed integrated circuits, it has not received as much attention as some of the other 
discontinuities, e.g., open-end terminations, bends, and junctions. This is due mainly to 
the nonplanar and complex 3-D geometry of the via, which has often been simplified in the 
works published previously [l]-[3]. For instance, a via penetrating through a single 
reference (ground) plane with two wire traces has been considered in [1], and without any 
traces in [2], while a via above a reference plane with two wire traces but without a 
through-hole reference plane between these traces has been investigated in [3], A novel 
equivalent network model, which accounts for the frequency dependence, has been 
proposed in [4] and has also been applied to the problem of coupling between two adjacent 
vias in [5]. In [1] and [3], an integral equation has been formulated in terms of the excess 
charge distribution to compute the equivalent (excess) capacitance. In this chapter, this 
excess charge formulation is further generalized for vias with more complex geometries 
than has been analyzed hitherto and is applied in conjunction with the closed-form Green’s 
function to analyze vias embedded in multilayered dielectric media. 

The general description of a via structure is discussed in Section 5.2. An integral 
equation is formulated using the closed-form Green’s functions in Section 5.3 to determine 
the excess charge distribution of a via, and the method of moments (MoM) is subsequently 
employed to determine the unknown charge distribution. A detailed discussion of the 
closed-form Green’s functions and the corresponding expressions can be found in Chapter 
2. Unlike the strip discontinuities considered in the previous chapter, the conductor traces 
associated with a via can be located at different vertical locations and a via often goes 
through several ground planes; hence, the formulation of an integral equation is slightly 
more complicated than the one given in the previous chapter. Although two different 
integral equations are used for a via with and without a through-hole ground plane in the 
other published methods [1], [3], only one integral equation is formulated to handle both 
cases in this chapter. In Section 5.4, several numerical examples are presented to verify the 
proposed method. 
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5.2 General Statement of the Problem 

To illustrate the geometries of vias considered in this chapter, a via comprised of 
three traces and one reference ground plane is shown in Fig. 5.1. In general, a via can 
pass through N g reference (ground) planes and N ( traces, and N p pads can be attached to 
the via where N g , N t , and N p are all arbitrary. The vias are embedded in a multilayered 
medium consisting of Nj (arbitrary) dielectric layers, which can be backed by two optional 
reference planes as shown in Fig. 5.2. To distinguish between these optional reference 
planes and those associated with the via, we will reserve the term “reference plane” to 
designate an optional top or bottom reference plane, and use the term “reference conductor” 
to denote other reference planes. It is evident that the reference conductors must have 
perforations to avoid any contact with the vias; however, the two optional reference planes 
are assumed to be solid. To simplify the numerical computation, we assume that all of the 
conductors are infinitely thin except a via hole, which is assumed to be filled with a 
conducting material, and that the shapes of all the pads and perforations in the reference 
conductors are circular. 

Two quasi-static equivalent circuit representations of the via shown in Fig. 5.1 are 
given in Figs. 5.3(a) and 5.3(b). These representations assume that there is at least one 
reference conductor or one reference plane. The two circuits are equivalent in the sense that 
one can be obtained from the other; however, the equivalent circuit shown in Fig. 5.3(b) is 
preferable because equivalent inductances, in general, are computed from terminal (trace) to 
terminal (trace), and these inductances directly correspond to the ones shown in Fig. 
5.3(b) 1 . This chapter will only address the problem of computing the total equivalent 
capacitance C e . The method to compute the equivalent inductance of a via can be found in 
[ 6 ]. 

5.3 Formulation of an Integral Equation 

An impressed potential on the conductors results in free charge accumulation on the 
surfaces of conductors, and the electrostatic potential <t>(r) at any point except inside the 
conductors is then related to this surface charge density q{r) via the following integral 
equation: 


'Note that for two-terminal (two-trace) case this is irrelevant since L e j and L e 2 are equal to L e j2> 2 . 
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<t>(r) = J G 3D {r\r' )q(J )dt> = (c 3D ,q) (5.1) 

n 


where £2 denotes the surfaces of all conductors, including the via and reference 
conductors. G (rlr 1 ) is the 3-D closed-form Green’s function for a multilayered 
medium, and it accounts for polarization charges on dielectric interfaces and free charges on 
the surfaces of reference planes. The integration is symbolically written as (•,•) to simplify 
the notation. Next, we represent the charge density q(r) as follows: 



(5.2) 


where q v (r) is the charge density on the surface of a via hole, and q l p (r), q l t (r ) , and q‘ g (r) 
are the charge densities on the ith pad, trace, and reference conductor, respectively. 
Equation (5.1) can be rewritten as 




N, 


N. 


«r) = (c 3 °« v ) + 2(G 3C ,«') + X(c 3D ,,;) + 2;(c 3D ,,') (5.3) 

1=1 1=1 1=1 

Next, the charge density q\{r) is decomposed into the uniform charge density 
qUnif,i( r ) and excess charge density qf xcess,l (r): 


q i t (r) = qr f ' i (r) + qf xcess ' i (r) 


(5.4) 


Here, q^'\r) is the uniform charge density on the ith trace under the assumptions that it 
is infinite in both directions, no other traces are present, and the reference conductors have 
no perforations; this charge density is computed by solving an appropriate 2-D problem. 
Since the reference conductors become uniform planes without any perforations for this 
2-D problem, the potential distribution on the region above the reference conductor is not 
affected by the region below it and vice versa. To solve for q^ if,i {r), it is then expedient 
to introduce a new medium surrounding the ith trace. As a consequence, the medium 
employed in the 2-D problem is generally different from that of the 3-D via problem, and 
could, in fact, be different for each trace. Once the appropriate medium has been chosen, 
qunif,i( r ) can be obtained by using the method described in Section 3.2. The resulting 
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qf n ‘f’i(r) yields the capacitance per unit length for the ith transmission line in the 
equivalent circuit representation shown in Fig. 5.3. 


In the process of determining q^ n ^'\r), the 2-D closed-form Green’s function 
G 2D 'i(p\p 0 ) is used to formulate an integral equation for a 2-D problem in Section 2.2. 
However, in the integral equation (5.3) for computing the equivalent capacitance problem, 
the uniform charge density q^'\r) resides on the ith trace, which is only a semi-infinite 
line. It is therefore necessary to employ G semi ' 1 (r\r 0 ,%) to compute the potential due to 
qUnif,i( r ). Using (5.4), (5.3) can be rewritten as 


N, 




i=i 


N„ 


V, 


N„ 


(G 3D . 9v ) + X(c 3D .?^ + £(c 3£> .9r w ')+X(G 3D .4) (5.5) 

1=1 i = 1 i=l 


Next, if we set the via potential to be 4> 0 with respect to the reference conductors and 
planes, <t>(r) becomes p 0 on the surfaces of the via hole, pads, and traces, and is equal to 0 
on the surfaces of reference conductors. Hence, once q^ ni ^'\r) has been determined, all 
of the quantities associated with the left-hand side of the above equation can be considered 
to be known at the surface of the conductors. 


The integral equation (5.5) can now be solved by using the method of moments. In 
this method, the surfaces of the conductors, with the exclusion of the reference planes, are 
first discretized with polygonal patches. Next, the unknown charge distributions, q v (r), 
q l p (r), qf xcess ' l (r), and q l g (r), are expanded with pulse-type basis functions over these 
patches. For instance, we wrote q v (r) as 

^v(r) = ^q v jPj(r) (5.6) 

;=i 


where Pj(r) is 1 if r is in the jth patch used to discretize the via hole and, 0, otherwise. 
Evaluating (5.5) at the centers of all the patches results in a system of linear equations, 
which yields the unknown charge distribution similar to (4.6) in Chapter 4. The various 
integrations appearing in (5.5) can be evaluated analytically for pulse-type basis functions 
using closed-form formulas given in Sections 2.3 and 2.4. 
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Once the unknown charge distributions have been determined, the equivalent 
(excess) capacitance C e can be obtained by using the following expression which involves 
the integrals of these charge distributions: 

f S N > 

C 'to = \q' p (ndr‘+Y, )ti“‘ SS ‘<r’)dr- (5.7) 

Q, i=] O p i 1=1 Q t . 

where Q v is the surface of a via hole, and Q pi and Q p i are the surfaces of the rth pad 
and trace, respectively. 

5.4 Numerical Examples 

First, two numerical examples are considered to illustrate the application of the 
method presented above to the computation of the equivalent capacitances of two via 
structures, one with a reference plane and the other with a reference conductor. The 
detailed geometries of the two via structures are shown in Figs. 5.4 and 5.5. The 
computed excess capacitances for Fig. 5.4 with 1) ej = e 2 = e 0 and 2) e } = 4e 0 and e 2 = £ ; 

for Fig. 5.5 with 3) £; = e 2 = 4e 0 and 4) £; = 4.5£ G and e 2 = 5.4e 0 are listed in Table 5. 1 

along with data obtained from [1] and [3]. In [1] and [3], the strips were replaced by the 
equivalent wires of radii which are one-fourth of the widths of the strips. In our 
computation, the lengths of all traces have been truncated to 2.5h, whereas the width of the 
reference conductor has been truncated to 1 .51, with h and / being the height of a via hole 
and the length of the traces. The truncation of traces and reference conductors is valid since 
the excess charge distribution decays rapidly as we move away from the center of a via. A 
total of 263 and 687 unknowns were used for the vias shown in Figs. 5.4 and 5.5, 
respectively. As shown in Table 5. 1, the data for the via shown in Fig. 5.4 agree well with 
the published results. However, the data for the via shown in Fig. 5.5 are considerably 
different from the results reported elsewhere. Unfortunately, no experimental result for 
this structure is available to establish the relative accuracy of these results associated with 
Fig. 5.5. 
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Reference Plane 


Figure 5.4. A two-trace via with a reference plane. All dimensions are in millimeters. 



Figure 5.5. A two-trace via with a reference conductor. All dimensions are in millimeters 
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Table 5.1. Equivalent capacitances for vias shown in Figs. 5.4 and 5.5. Units are in 
picofarads. 



Fig. 5.4 


Fig. 5.5 


e l= e 2~ £ o 

£y— 4 e o , £2"“^o 

E] =£ 2 =4£ o 

£/=4.5 £ 0 , £ 2 =5.4£ 0 

Computed Result 

0.3841 

1.233 

9.952 

12.31 

Others 

0.3701 [3] 

1.28 [3] 

6.35 [1] 

7.85 [1] 


The equivalent capacitance of the via shown in Fig. 5.6 is obtained from the 
experiment based on Hewlett Packard Application Note #67 using a 500 ps rise-time input 
signal; the resulting capacitance value was 0.504 pF. To apply the presented method to this 
via, only that portion between the ground planes is modeled in the computation, and it is 
assumed the two ground planes are uniform. The total of 382 unknowns is used in the 
computation and the computed excess capacitance value was 0.395 pF. 

The comparisons of the CPU time at the Alpha workstation for various matrix 
solution techniques are given in Table 5.2. Again, the GCR routine performed best for all 
examples as in the previous chapter. 

5.5 Summary 

A method to compute the equivalent capacitance of a via, which is based on an 
integral equation formulated in terms of the excess charge formulation, has been presented 
in this chapter. The method is applicable to via geometries with or without through-hole 
reference conductors. The closed-form Green’s function was employed to circumvent the 
time-consuming evaluation of a nested infinite series, required in the evaluation in [3]. 

Additional computational savings can be achieved via the use of the Fast Multipole 
Method (FMM) [7], an algorithm to speed up the computation time in 3-D capacitance 
calculations, in conjunction with the closed-form Green’s function for a multilayered 
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Figure 5.6. A two-trace via used to obtain the experimental data. Units are in millimeters. 


Table 5.2. The CPU time used in the computations. 



GCR 

GMRES 

LU 

Time (sec) 

Iterations 

Time (sec) 

Iterations 

Time (sec) 

Fig. 5.4 

13.18 

28 

14.14 

52 

17.12 

Fig. 5.5 

81.22 

35 

120.08 

50 

161.93 

Fig. 5.6 

167.46 

12 

167.52 

17 

171.76 
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medium. In the following chapter, such an implementation of the multipole method with 

the closed-form Green’s function is demonstrated for the computation of the equivalent 

capacitance of a strip crossover. 
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CHAPTER 6 

COMPUTATION OF THE EQUIVALENT CAPACITANCE OF A STRIP 
CROSSOVER IN A MULTILAYERED MEDIUM 

6.1 Introduction 

In Chapters 4 and 5, the equivalent capacitances of various junction discontinuities 
and a via have been considered based on the quasi-static approximation. In this chapter, 
the equivalent capacitance of yet another common discontinuity, a strip crossover, is 
considered. Like the methods presented in the two previous chapters, the method 
discussed in this chapter is also based on the quasi-static approximation and utilizes a 
closed-form Green's function. 

To compute the equivalent capacitance for orthogonally crossing strips without a 
top ground plane, an iterative spectral-domain approach is used in [1], a spatial Green's 
function with the infinite series expansion is used in [2], and a spatial Green’s function 
based on the complex images is used in [3]. A crossover embedded in lossy multilayered 
media is also considered in [4], The full- wave analysis of a strip crossover appears in [5]- 
[7]. In this chapter, the closed-form Green’s function is applied to a crossover of an 
arbitrary crossing angle and an additional top ground plane. 

The Fast Multipole Method (FMM) [8]-[12] is a class of algorithms used for the 
rapid evaluation of potentials in large systems The techniques in [8]-[l 1] have been applied 
to compute the capacitance of conducting structures in a homogeneous medium [10] and in 
the presence of finite-sized dielectrics [14]. Anderson's fast-multipole method [12], which 
is equivalent in complexity and accuracy to the methods described in [8]-[ll], has been 
applied in conjunction with the closed-form Green’s function for computing the capacitance 
matrices of conducting structures that reside in a layered dielectric medium [15]. In this 
chapter, the possible usage of the FMM to accelerate the MoM computation of the 
equivalent capacitance calculation is demonstrated. Numerical experiments show that a 
substantial amount of computer time is saved without a loss of accuracy. 
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6.2 General Statement of the Problem 

Consider two strips crossing at an angle a as shown in Fig. 6.1(a). The bottom 
and top strips are denoted by Conductors 1 and 2, respectively. Although this chapter will 
focus on strips, the proposed method is applicable to conductors of arbitrary cross section. 
The strips are embedded in a multilayered dielectric medium depicted in Fig. 6.1(b). An 
arbitrary number Nj of dielectric layers are located on top of a ground plane, and the 
multilayered dielectric medium is terminated by an optional ground plane on the other side. 
The dielectric medium is assumed to be uniform and of infinite extent along the x- and y- 
directions. 

When compared to bends or junctions, the discontinuity effect associated with a 
strip crossover is less localized due to the absence of the nonuniformity of conductor 
geometries. As a consequence, an equivalent circuit with distributed parameters may be 
needed to accurately model the crossover. Although the required parameters can be 
obtained by using the proposed method, a circuit with distributed parameters is not 
convenient to simulate in practice. Only the simple lumped equivalent circuit shown in Fig. 
6.2 will be used in this chapter. Unlike the discontinuities considered in the previous 
chapters, a crossover consists of two electrically isolated conductors; hence, the equivalent 
circuit representation of a crossover consists of three capacitances, and a coupled integral 
equation has to be solved to determine these capacitances. In general, an equivalent 
inductance should be incorporated into the equivalent circuit to model current coupling for 
nonorthogonal crossovers; however, this chapter only concentrates on the computation of 
an equivalent capacitance. 

The remainder of this chapter is organized as follows. In Section 6.3, a coupled 
integral equation in terms of excess charge densities is formulated using a closed-form 
static Green's function described in Chapter 2. Then, the Fast Multipole Method is 
introduced in Section 6.4, and it is demonstrated how it can be used to accelerate the MoM 
solution of the integral equation. Finally, various numerical examples are presented in 
Section 6.5. 
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Figure 6.2. Equivalent circuit representation of a strip crossover. 

6.3 Formulation of an Integral Equation 

Let us first assume that the strip crossover depicted in Fig. 6. 1 (a) is embedded in a 
multilayered dielectric medium with both top and bottom ground planes. An impressed 
electrostatic potential on conductors results in charge accumulation on the surfaces of 
conductors. The electrostatic potential *f^r) at any point is related to this accumulated 
surface charge density q(r) through the following integral equation: 

nr) = j 0 G 3D (rSMW =(G 3D ,?) (6.1) 

where G 3D is the 3-D Green's function for the multilayered medium, and Q denotes the 
surface of the conductors. The integration is again symbolically written as {•,•) to simplify 
the notation. 

As noted in Chapter 2, the spectral-domain Green's function for a layered medium 
with both top and bottom ground planes cannot be approximated with a finite number of 
exponential functions owing to the existence of a pole; therefore, both 2-D and 3-D Green’s 
functions in the space domain cannot be represented by a finite number of weighted 
images. Although the exact closed-form expression (without weighted images) can still be 
found for the 2-D space-domain Green's function (see Appendix B), such a closed-form 
expression cannot be found for the 3-D space-domain Green's function; the expression of 
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the 3-D Green's function consists of slowly converging infinite series due to the infinite 
number of image charges. To avoid the use of this 3-D Green's function, the top ground 
plane is truncated and modeled as an additional conductor; hence, G 30 used in this chapter 
refers to the 3-D Green's function for a multilayered dielectric medium with only the 
bottom ground plane. This approach results in a large number of unknowns due to 
modeling the top ground plane. It may therefore be computationally inefficient to use this 
approach as compared to using the Green’s function with an infinite series expansion. 
However, when a multipole-accelerated algorithm is applied, it is expected that the 
proposed approach will result in better computational efficiency in terms of memory usage 
and CPU time than that due to the use of the Green’s function with an infinite series 
expansion. This is primary because of the small number of images required in the closed- 
form Green’s function. 

Now let us represent the surface charge density q(r) as follows: 

q( r ) = Qc \ ( r ) + <7c2( r ) + <lTplane(. r ) ($.2) 

where q c \{r), q c 2 (r), and q Tp iane( r ) 316 surface charge densities on Conductor 1, 
Conductor 2, and the top ground plane, respectively. The Green's function G 3D account 
for the charge density on the bottom ground plane and the polarization charge densities on 
the dielectric interfaces. By exciting Conductor 1 to a voltage <j> 0 with respect to Conductor 
2 and ground planes, (6.1) becomes 

<t> 0 =(G 3D ,ql l ) + (G 3D ,ql 2 ) + (G 3D ,qT p iane) for r on Conductor 1 (6.3a) 

0 = (g 3D , q\\ ^ + (G 3D ,q X c2 j + (g 3D >qTpUme) 

for r on Conductor 2 and top plane (6.3b) 

The superscript 1 in q\ x , q X c2 , qjplane denotes the charge densities due to the excitation of 
<j) 0 on Conductor 1 with respect to Conductor 2 and the ground planes. Now let 

i'c2 


(6.4b) 
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where q^\ if) is the uniform charge density obtained by solving the 2-D problem for 
isolated Conductor 1, i.e., without Conductor 2 but with the top ground plane. This 
problem is solved using the 2-D Green's function G 20 with the excitation voltage set to <p 0 . 
A detailed discussion of this procedure is given in Section 3.2. Equation (6.3) can then be 
written as 


to = <G 2D , 9 ;r r ) + {c 3 Vl“"“) + {o }D .q' c f c,ss ) + (c 3D ,q)f‘ p Ze) 


0 = ( G 2D ,q ' c ]“ n ‘ f ) + (0 3D ,q' c f““) + (G 3D ,«;f““) + (o W ,q 


excess ^ 


(6.5a) 

(6.5b) 


Here, is the 2-D Green's function for a layered medium with both top and bottom 
ground planes, and G 30 is the 3-D Green's function for a layered medium with only a 
bottom ground plane. 1 In (6.3) and (6.5), qj-piane and q represent two different 
charge distributions on the top ground plane: qj-piane is the induced total charge distribution 
due to the total charge distributions on the strips, whereas q^^e * s Educed charge 
distribution due only to the excess charge distributions on the strips. The induced charge 
distribution on the top ground plane from the uniform charge distribution on the strips is 
embedded in G 20 in (6.5). Since both the excess charge distributions on the strips and the 
induced partial charge distribution on the top ground plane from these excess charge 
distributions are all localized around the cross-over region, the strips and a top ground 
plane can be safely truncated when modeling the excess charge density q^u^e ■ 

Noting that (g 2D ,ql'“ ni ^\ is equal to (j) 0 on Conductors 1 and 0 on the top ground 
plane, (6.5) becomes 

(\_lr3D l,excess\,/ / ^3D l,excess\ , //-3D A,ex cess\ 

°-{ G ’<7cl ) + \ G ) + \ G biplane ) 

for r on Conductor 1 and top plane (6.6a) 

-(c 2 D ,q'r f ) = (o iD ,,;r”) + ( g 30 ,,;-^) + (g 30 ,,^) 

for r on Conductor 2 (6.6b) 


'When the top ground plane is not present, both G 2D and G 3D arc the Green's functions for the 
same medium, and the exponential approximation can be performed only once in the spectral domain to 
obain expressions for both G 2D and G- . 
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Following a similar procedure, but exciting Conductor 2 to a potential <p 0 with 
respect to the other conductors instead of Conductor 1, leads to 

/ r^2D 2,unif\ / y~,3D 2,excess\ , / ^D 2 ,excess\ , I f3D 2,excess\ 

~\ G '<* c2 / = \ G '4 cl / + \ G '4c 2 l + \ ij Tplane j 

for r on Conductor 1 (6.7a) 

n I r-^D 2,excess\ .! r 3D 2 ,excess\ , I r 3D 2,excess\ 

0 = \ G .?cl / + \ G 'Qc2 / + \ G Tplane j 

for r on Conductor 2 and top plane (6.7b) 


The superscript 2 in q 2 { excess , q^™ 6 ™ > denotes the excess charge densities as a 

result of the excitation of Conductor 2. 


The method of collocation, which is based on the pulse-type basis functions and the 
delta testing function, is used to solve the above coupled integral equations, (6.6) and 
(6.7), and the resulting linear system of equations is 


0 V 21 

V12 0 

0 0 


l,unif 

q l 

0 


where 


0 

2,unif 

<*2 



Mj 2 

m 13 * 


*.2, excess’ 


M 2 i 

m 22 

M23 

..1, excess 
q 2 

2, excess 

<*2 

(6.8) 

M 3 i 

M 32 

M 33_ 

1, excess 
“Tplane 

2, excess 

“Tplane 




^ds' 



(6.9a) 




(6.9b) 


-|7 

, excess 1 , excess 

,2 


i, excess _ f i, excess A>e 

q j ~L q j' x ,qj ' 7 

i, excess _ \ Lexcess i y excess i, excess 1 

^Tplane ^Tplane.V^Tplane^'^^Tplane.Njp^ J 

T 

„i,unif _ \ -2D 2D 2D "1 
Qi — l^i.l ’9i,2 J 


(6.9c) 

(6.9d) 


(6.9e) 
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Here, S q and C q denote the source patch and line segment for 3-D and 2-D problems, 
respectively. Nj is the total number of patches used to represent q l j excess (r), and N?° is 
the total number of line segments used to represent q£ ni ^ (r). Nppiane is the total number 
of patches to represent g f j^^ s (r) . It should be emphasized that all quantities appearing in 
the left-hand side of (6.8) are known assuming that is precomputed by solving the 
2-D problem for isolated Conductor i. 

The equation corresponding to (6.8) for a case without the top ground plane can be 
simply expressed as 


' 0 

V 2 f 

■ l,unif 

0 



m 12 

1, excess 

_ 2, excess" 
Ml 

. V 12 

0 

0 

2, unit 

q 2 J 


M 21 

m 22 . 

_1, excess 
_ q 2 

^2, excess 
q 2 


Once (6.8) or (6.9) is solved, the equivalent capacitance can be obtained by 


excess 

H 

^\, excess 
L 2 

r 2 , excess 1 f f\\ y ex cess 

C 1 _ 

^2, excess sA ^excess 

c 2 j ifj 2 

s\2, excess' 

s{L,excess 

\l2 

(6.11a) 

c m 

_ ^1, excess ^-.2, excess 2 

- -t- 2 - -Cj 

(6.11b) 


C? _ Qiy excess _ Qm 


(6.11c) 


where 

«J 

Q? xcess = X qi i"“ eSS ' Area J- l ‘ ( 6.1 id) 

k=] 

and Areajjc is the area of the kth patch in conductor j. In the above formulas, it is assumed 
that (6.8) or (6.10) is solved with <j) 0 set to 1. 

It should be noted that the difference between the integral equations for the general 
3-D capacitance computation and the equivalent excess capacitance lies only in the left side 
of (6.8) or (6.10), i.e., in the excitation vectors. This fact allows the Fast Multipole 


2 The later equality is due to the reciprocity; however, C'“"“ and C 2 '""" , in general, are slightly 
different due to the numerical side effect, and the average scheme should be used to compute C”. 
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Method (FMM) in [15] to be employed in solving (6.8) and (6.10), as shown in the 
following section. 

6.4 The Fast Multipole Method (FMM) 


The solution of (6.8) can be obtained through solutions of the equation having the 

form 


V = MI 


(6.12a) 


where 


V = 


0 

v u q*' 

0 


(6.12b) 


and 


1 = 


q:“"“ 

q 

qST 


(6.12c) 


or 


and 


V* < 11 °°“ 



(6.1 2d) 


1 = 


n 2,«cta 

qr 

qiocw 


(6.12e) 


and 



114 



Equations (6.8) and (6.10) can be solved either directly by LU decomposition, or 
by iterative techniques such as the Generalized Conjugate Gradient (GCR) or the 
Generalized Minimum Residual (GMRES). If the MoM matrix, represented by M in 
(6.12a), has dimensions NxN, then a direct method, for example LU decomposition, 
would require 0(N 3 ) operations. If an iterative technique is used to solve the MoM 
system, the cost would be 0(pN 2 ), where p is the number of iterations since a matrix- 
vector product is required in each iteration. This product has the form MI^ where l q is a 
vector of trial charge coefficients in the gth iteration. The FMM can be invoked to compute 
this product in an extremely efficient manner, requiring only 0{M) operations per iteration, 
as demonstrated in [12]. This results in an overall complexity for the FMM-based iterative 
technique of 0(pM), where p is the number of iterations [13], [14). 


The approach followed here, and explained in detail in [15], is to use Anderson's 
FMM technique, which is based on the use of Poisson's formula to represent the solution 
of Laplace's equation. The images resulting from the use of the closed-form approximation 
to the 3-D Green's function G 3D has to be accounted for in the FMM. To facilitate the 
required computation, the entire multiconductor structure and all image patches are enclosed 
in a large cubical volume (called the parent cube), which is recursively partitioned into eight 
smaller cubes (termed the child cubes). At the nth level, the problem space has been 
hierarchically partitioned into 2 3n cubes. For each cube at every level, an "outer sphere" 
approximation is constructed [12], [15]. This approximation permits the computation of 
the potential at any point outside the sphere, given the potential on its surface due to all 
enclosed sources. These approximations are constructed in a hierarchical manner to 
efficiently aggregate source points. Observation points are aggregated in a manner similar 
to that employed to aggregate source points, by using "inner sphere" approximations that 
facilitate the computation of the potential at any point inside a sphere, given the potential on 
its surface from all sources outside it. The surfaces of the conductors are divided into 
genuine patches, while the images produced by the Green's function lead to the presence of 
image patches. The observation points correspond only to the centers of the genuine 
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patches, whereas both genuine and image patches comprise the sources: the potential 
boundary conditions are enforced only on the conductor surfaces. The overall FMM-based 
technique for evaluating the matrix-vector product requires an initial potential evaluation 
due to every source, construction of outer and inner sphere approximations, and a final 
potential at evaluation at every observation point. Note that given a charge distribution, the 
FMM obtains the potential values without ever generating the matrix M. 

At present, our FMM implementation works for the case of a single dielectric layer 
backed by a ground plane, with free space on top. Finite ground planes can be modeled. 
The extension to an arbitrary number of layers is relatively straightforward and is currently 
under implementation. 

6.5 Numerical Examples 

In this section, the usefulness of the above algorithms in the evaluation of the 
excess capacitance for various representative geometries is demonstrated. It is assumed 
that all strips considered in this section have negligible thicknesses. Since the excess 
charges are localized around the discontinuity, the lengths of all strips are truncated to the 
finite length of / using the following formula: 

/ = 20 (/*2 -/*!)— (6.13) 

a 

where h i and /12 are the heights of Conductors 1 and 2 from the bottom ground plane, and 
oc is the crossing angle in degrees. The length of a top ground plane lg, if it exists, is taken 
to be 21. Nonuniform meshing is utilized to descretize both conductors and the top ground 
plane. Unless otherwise specified, the total numbers of basis functions for 2-D and 3-D 
problems are 16 and 160 for each strip, respectively, and 625 basis functions are used to 
model the top ground plane. A maximum of five exponential functions were used to 
approximate each coefficient function of the Green's function, and the Generalized 
Conjugate Gradient (GCR) is used to solve the linear system of equations. The CPU time 
is measured on the Alpha workstation. 

As a first example, consider the strip geometry depicted in Fig. 6.3. The width of both 
strips is 0. 16 mm. The equivalent capacitance is computed over a wide range of crossing 
angles, and the results are shown in Fig. 6.4. The excess charge distribution is plotted in 
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Conductor 2 Free Space 



Bottom Ground Plane 

Figure 6.3. Example 1: a microstrip crossover. 


Equivalent Capacitance vs. Crossing Angle 



= 6 mm 


Figure 6.4. Variation of the equivalent capacitance as a function of the crossing angle a. 
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Fig. 6.5. Table 6.1 shows the comparison data for the orthogonal crossing case. The 
difference between data in the first and second columns is due to the difference in meshing 
size and the order of approximation used to obtain the closed-form Green's function. A 
total of 1 1.36 sec of CPU time is used to compute the orthogonal case. 


Next, consider the crossover structure with a top ground plane shown in Fig. 6.6. 
The width of the strips is again 1 mm. Figure 6.7 shows the variation of the values of the 
equivalent capacitance over the height of a top ground plane hj. The comparison of the 
result for /ij = 50 with that for the case without a top ground plane is shown in Table 6.2. 
A total of 55. 1 1 sec of CPU time is used for the case with a top ground plane. 


Equation (6.1) is also solved directly for the total charge distribution 3 instead of 
formulating the integral equation in terms of the excess charge distribution, which is the 
conventional way to find 3-D capacitance matrices. Then, the following formula is used to 
compute the excess (equivalent) capacitance matrix: 


C l, excess r~l2 ,excess 

1 C 1 

fl, excess ^2, excess 

c 2 c 2 


q\, total _ I * q1, uniform 

Q l, total 
2 


n 2, total 

^1 

r,2,total i * n 2,uniform 

-tty 


(6.14) 


where Q\ ,total is the sum of the total charge distribution on Conductor j with the excitation 
on Conductor i, and Q^ uni f orm i s the sum of the uniform charge distribution on Conductor 
i. As mentioned in Section 4. 1 , the major drawback of this approach is that since the 
excess charge due to the crossover is much smaller than the total charge, the final accuracy 
in terms of the excess charge is much worse than the accuracy obtained for solving the total 
charge. Table 6.2 demonstrates this fact. 

For the final example, the crossover shown in Fig. 6.3 with hj = 2.5 mm, /12 = 5 
mm, w] = W 2 = 1 mm and e r = 5 is considered to demonstrate the application of the Fast 
Multipole Method (FMM) described in [9], [12]. The FMM is used with the Generalized 
Minimum Residual (GMRES) routine, as in [12]. The comparison, in terms of the CPU 
time, with the LU factorization, the Generalized Conjugate Residue (GCR) routine, and the 


term total charge denotes the sum of uniform and excess charges. 



Excess Charge Density (C/(mm x mm)) Excess Charge Density (C/(mm x mm)) 
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Figure 6.5. Plots of the excessive charge distributions for (a) qff, (b) gfi> ( c ) 20(1 

(d) <722 f° r Example 1 for the orthogonal crossing case. 
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Table 6. 1 . Comparison data for the orthogonal case of Example 1 . 



Computation 

Data from [2] 4 

Data from [3] 

c m 

64.81 fF 

65.12 fF 

65.16 fF 

cf 

-51.30 fF 

-51.42 fF 

-51.58 fF 

ci 

-55.10 fF 

-55.08 fF 

-54.92 fF 


Top Ground Plane 



Bottom Ground Plane 

Figure 6.6. Example 2: a strip crossover with a top ground plane (orthogonal crossing 
case). 


‘‘The entries in Table 2 in [2] contained errors, and the corrected data were obtained from the 
authours and given in Table 6. 1 . 



121 


<D 

O 

a 

B 

'§ 

8 * 

u 


The Equivalent Capacitance vs. the Height of 
a Top Ground Plane 



Figure 6.7. Variation of the equivalent capacitance as a function of the height of a top 
ground plane. 


Table 6.2. Comparison data for the orthogonal case of Example 1 . 



Without Top Ground 

With Top Ground (/ij = 50) 


The Excess Charge 
Formulation 

The Total Charge 
Formulation 

The Excess Charge 
Formulation 

c m 

332.75 fF 

309.87 fF 

329.46 fF 

^9 

-291.23 fF 

-139.40 fF 

-286.20 fF 


-298.24 fF 

37.11 fF 

-296.65 fF 
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Computation Time Comparison 



Figure 6.8. Comparison of CPU time used for the LU factorization, GCR, GMRES, and 
FMM. 


Generalized Minimum Residual (GMRES) routine, is given in Fig. 6.8. The result from 
LU factorization was C m = 518.06 fF, Cf = -415.88 fF, and C| = -460.98 fF. The 
FMM gave C m = 517.0 fF, Cf = -414.4 fF, and C| = -460.1 fF. Although this particular 
numerical example can be solved without the FMM, in general, the number of unknowns 
can be quite large when modeling the thicknesses of the strips and the top ground plane, 
and the usage of FMM for such cases will significantly improve the computation speed as 
indicated by Fig. 6.8. 

6.6 Summary 

The computation of equivalent capacitances of a strip crossover is considered in this 
chapter. Strips crossing at an arbitrary angle with both top and bottom ground planes are 
considered for the first time in this chapter. The presented method is based on a static 
integral equation in conjunction with a closed-form Green's function, and does not involve 
numerical integrations or infinite summations in the evaluation of the MoM matrix. The 
application of the Fast Multipole Method to accelerate the computation is also discussed. 
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CHAPTER 7 

SUMMARY AND FUTURE WORK 

In this study, computationally efficient and accurate methods to compute the 
transmission line parameters of interconnections and the equivalent capacitances of 
discontinuities associated with interconnections were presented based on the quasi-static 
approximation. An extensive amount of research has already been performed on this 
subject and various numerical techniques have been developed thus far, as referenced 
throughout the chapters; however, most of these techniques are applicable to only a certain 
class of problems or are computationally intensive. For instance, some methods assume 
conductors to be of zero thickness and are limited by the number of dielectric layers (often 
two layers), whereas some methods require an evaluation of nested infinite summations, a 
numerical integration, or modeling of additional unknowns over dielectric interfaces. In 
contrast, the presented methods are applicable to multiple conductors embedded in 
multilayered dielectric media, and the cross sections of conductors can be arbitrary 
although, in the cases analyzed, we have used infinitely thin strips for modeling of 
discontinuities to simplify the analysis. Moreover, these methods are numerically efficient 
as they do not involve any of the computationally expensive operations described before. 

The computational efficiency of the presented methods is mainly due to the 
introduction of the closed-form Green’s function for a multilayered dielectric medium. 
Although the concept of the closed-form Green’s function first originated in full-wave 
analysis and had already been applied to a static case using complex images, the closed- 
form Green’s function based on real images, which is computationally more efficient than 
the method based on complex images, was first developed in this study. In fact, the 
closed-form Green’s function (based on complex images) for static analysis was published 
at an early stage of this study, and a similar concept using real images was independently 
inspired by this author. The complete list of expressions for the closed-form Green’s 
functions for a point charge, a uniform line charge, and a semi-infinite uniform line charge, 
all embedded in a multilayered dielectric medium, was presented in Chapter 2 with the 
associated closed-form integration formulas for the collocation method. 

The closed-form Green’s functions derived in Chapter 2 were then applied to solve 
various electrostatic problems. In Section 3.2, this Green’s function was used to solve 2-D 
electrostatic problems to compute the capacitance matrix of a multiconductor system. In 
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Chapters 4, 5, and 6, it was used to formulate integral equations in terms of the excess 
charge distributions to compute the equivalent capacitances of various interconnection 
discontinuities, such as open ends, bends, various junctions, vias, and crossovers. 

The electrostatic solution obtained in Section 3.2 was further used to compute the 
conductance matrix of a multiconductor system in Section 3.4. The normal component of 
the electric field at the surfaces of conductors was determined from the charge distribution 
obtained from Section 3.2; then, it was further related to the shunt current density at the 
surfaces of conductors and used in the computation of the conductance matrix. Hence, the 
presented method for the computation of the conductance matrix is perturbational in the 
sense that the normal component of the electric field at the surfaces of conductors of a 
lossless system was used to determine the shunt current density. 

To compute the inductance and resistance matrices, an analogy between electrostatic 
and magnetostatic problems was presented in Section 3.3. Then, an equivalent electrostatic 
problem was solved instead of a magnetostatic problem, and the inductance matrix was 
obtained from the capacitance matrix of the equivalent electrostatic problem in Section 3.3. 
The conduction current density on the surfaces of conductors was also obtained from the 
surface charge density of the equivalent electrostatic problerii, and it was used to determine 
the resistance matrix in Section 3.5. To include losses due to imperfectly conducting 
ground planes, coupled integral equations were formulated in Section 3.5 to obtain the 
current densities on the surfaces of the ground planes. 

Throughout the chapters, the results obtained from the presented methods were 
compared with other published results including some experimental results and, in general, 
a good match was obtained. 

The Fast Multipole Method (FMM) is a recently developed algorithm to accelerate 
the evaluation of potentials in large system. Some authors applied this multipole algorithm 
to a multilayered case by modeling the polarization charge on dielectric interfaces; however, 
a more efficient implementation of the multipole algorithm to a multilayered case may be the 
use of the closed-form Green’s function. At the present time, the multipole algorithm has 
been applied to the closed-form Green’s function for only a stratified medium where all 
conductors are located in the same dielectric layer. In this thesis, this method was 
considered in the computation of the equivalent circuit of a crossover to obtain further 
savings in the computation time and memory; a similar approach can also be applied to 
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computation of the other discontinuities. The implementation of the multipole algorithm to 
a more general closed-form Green’s function, in which conductors can be located in any 
layer, would be a good subject for future work. 

Although only the capacitive nature of the discontinuity is considered in this thesis, 
in general, the effect of the discontinuity will be far more complex and full-wave analysis 
must be employed to circuits in the high frequency regions. The author is currently 
investigating the application of the Finite-Difference Time-Domain (FDTD) method to 
generate the equivalent circuits of discontinuities to account for the frequency-dependent 
nature of these nonuniform structures. 
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APPENDIX A 

THE REAL-VALUED EXPONENTIAL APPROXIMATION 
BASED ON THE RELAXATION OF CURVE FITTING 

First, we will assume that a function y(x) to be approximated is real valued and 
nonoscillatory and, furthermore, that its asymptotic value is zero. The latter assumption 
can easily be satisfied if the function is limited at infinity. Our goal is to find the right-hand 
side of 

N N 

y(x) = f a (x) = 2>) = ^qe X ‘ x (A.l) 

i=l i=l 

Let us for a moment assume that each first-order function f^x) approximates the 
original function y(x) at some interval around one of the approximation points and is 
decreasing fast enough so that its value is negligibly small at the approximation points 
corresponding to larger values of the argument x. Then, we can safely determine one of 
the first-order functions, say f\(x), by neglecting contributions due to the other first-order 
functions, which are unknowns to be determined. The parameters of f x (x) can be easily 
obtained by curve fitting two values of y(x) for some large value of x. In a similar 
manner, we can find the parameters of the other first-order functions; however, this time 
we have to take into account contributions due to the previous ones which are already 
known. 


From the above argument, given IN approximation points, the equations used to 
determine the parameters for the ith first-order function are then written by 


^ _ 27tln(y,(x 2t _ 1 )/y t (x 2 ,)) 

x 2i - *2,-1 


(A. 2) 


q=e A ' X2 '-iy l (x 2f _i) 


(A. 3) 


where 

yi ( x j ) = yi-i (*y ) - fi -\ ) ; j = i, -,2N (a.4) 


In the above, yq (•*./) equal to y(xj). Let us now consider the case in which the 
value of a first-order function is not negligible at the other approximation points. In this 
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case, if we perform the above procedure, there will be some difference between the original 
and the approximated functions since we have ignored contributions due to some of the 
first-order functions. In such a case, to reduce this difference, we can iterate the above 
procedure including the contributions from all other first-order functions which were 
obtained from the previous iteration. Thus, for the kth iteration, (A.4) must be modified as 

= j = \. ;2N (A.5) 

/=1 l=i + 1 

It can easily be shown that if y(x) had N distinct eigenvalues, by iterating the above 
procedure, the approximate function will converge to the original function. However, in 
general, the approximate function will never be exact and the iteration must be stopped at 
some point at which the approximate function is optimal in the curve-fitting sense. Since 
our curve-fitting algorithm, most likely, has the largest relative error for values of x 
between X 2 N- 1 and X 2 N, one can check the approximated value with the exact value at any 
point in this interval for the convergence criterion. In some cases, if the desired accuracy 
can not be achieved, then one should increase the number of exponentials. This case can be 
determined by checking the difference between the computed parameters of f- k) and 

The described method allows one to find the parameters of an approximating 
function one-by-one directly without solving a system of nonlinear or linear equations and 
does not require an original function to be monotonic. Moreover, the limiting values of the 
approximated function match exactly with the original function. 

Finally, to demonstrate the method, the following testing function is approximated 
and the results are shown in Fig. A. 1 : 

y(x) = -^— + e x (A.6) 

As illustrated in Fig. A. 1(b), to find the exact location of the pole at -1, we need a 
large number of iterations; however, since our goal is to approximate the overall function, 
only a few iterations were needed to approximate the function. To locate the poles exactly, 
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one should consider other methods such as those based on pencil of functions [1] or the 
Prony approximation [2]. 
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APPENDIX B 

THE GREEN’S FUNCTION EXPRESSION FOR A STRIP LINE 


Solving Poisson's equation with separation of variables, the Green's function for 
the strip transmission line case with the distance h between two ground planes (see Fig. 
B.l) can be written as follows (see [1]): 


(B.i) 

ne*-*n v h ) V h ) 

n = 1 


The above expression for the Green's function can be easily integrated analytically over the 
pulse basis function. This summation is quickly converging due to the decaying 
exponential factors, and only a few first terms are needed; for instance, four digits of 
accuracy can be obtained using less than five terms for \x- jc’|//i> 0.5. However, when 
\x — x'\/h is small, the summation converges rather slowly, and the use of this Green's 
function expression should be avoided. 


Fortunately, the above series can be summed up using the following formula: 



(B.2) 


and the resulting closed-form expression is given by [(1.462), 2] 


1 In 
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sinh 

k(x — x') 
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+ sin 2 

"rc(y + y )" 
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sinh 2 

'k(x-x')' 
. 2 h 

+sin 2 

"jt(y-y' )‘ 
- 2 h . 



(B.3) 


This expression for the Green's function can also be obtained using the conformal map ping 
and the method of images (see Ch. 10 in [3]). The major disadvantage of this expression 
compared to (B.l) and (2.24a) is that the closed-form integration over the pulse basis 
function is unappealing; hence, the numerical integration scheme is unavoidable. 
Therefore, this form of the Green's function must be only used when \x-x'\/h is small. 
As discussed in [4], to integrate (B.3) numerically, it is convenient to rewrite this 
expression by extracting the singularity as 
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Ground Plane 

— — — y=h 


' y=y' 


Ground Plane 


y= o 



Figure B. 1 . A strip transmission line. 


G(x, y\x' , y ) = ln[(x - x' ) 2 + (y - y’ ) 2 ] + g( x, y|x' , y' ) (B .4) 


where 


gU,;y|jc\y) = 



47te 


[(x-x' ) 2 + (y - y' ) 2 

sinh 2 

’7C(jc — JC* )~ 
2 h . 

+ sin 2 

'jt(y + y)T 
2 h 1 

sinh 2 

'k(x-x') 
2 h 

. 7 

+ sin 

'7t(y-y )' 
2h 



(B.5) 


Now the first term in (B.4) can be analytically integrated over a line segment using the 
formula given in [5], but the second term must be integrated numerically. 
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